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CHAPTER  ONE 


INTRODUCTION 

1.1  PROBLEM  DEFINITION 

An  underwater  acoustic  array  is  a  specific  spatial  arrangement  of 
a  set  of  transducer  elements.  It  can  be  used  as  either  a  transmitter 
or  a  receiver.  Each  of  these  modes  of  operation  places  certain  re¬ 
quirements  on  the  performance  of  the  array.  The  degree  to  which  these 
requirements  are  met  determines  the  effectiveness  and  usefulness  of 
the  array. 

The  directional  characteristics  of  such  an  array  are  controlled 
by  the  shading  values  at  each  position  in  the  array.  These  shading 
values,  or  weights,  are  developed  under  the  assumption  that  all  the 
elements  have  identical  performance  characteristics,  expressed  in 
terms  of  amplitude,  (or  gain),  and  phase  over  a  broad  frequency  range. 

The  variations  found  in  actual  elements  refute  this  assumption. 
The  introduction  of  element  errors  degrades  the  performance  character¬ 
istics  of  the  array.  Because  of  this,  the  side-lobe  levels  may  de¬ 
viate  significantly  from  the  design  specifications. 

The  amplitude  and  phase  characteristics  of  the  response  from  each 
element  are  also  influenced  by  the  element's  location  within  the 
array.  This  positional  dependence  results  from  the  acoustic  inter¬ 
action  with  all  the  other  elements  in  the  array. 


A  reasonably  effective  two  part  method  of  reducing  the  amount  of 
error  introduced  into  the  response  pattern  by  elements  with  different 
tolerances  was  developed  by  Kendig  [1].  Some  of  this  method's  draw¬ 
backs  have  been  solved  by  Schafer  {2).  In  addition,  Schafer  has  de¬ 
veloped  a  computer  version  of  this  method. 

Originally,  the  first  part  of  the  Kendig  method  used  indirect 
measurements  of  frequency  and  capacitance  to  estimate  the  amplitude 
and  phase  errors  of  each  element.  The  second  part  of  the  Kendig 
method  determined  the  placement  of  elements  in  the  array.  The  Kendig 
Scatter  Diagram  Method  involves  the  pairing  of  elements  in  a  manner 
that  partially  cancels  out  their  respective  errors.  With  the  develop¬ 
ment  of  an  in-air  current/velocity  measurement  technique,  Schafer  [2] 
was  able  to  use  a  lumped-parameter  equivalent  circuit  to  predict  more 
accurately  the  amplitude  and  phase  characteristics  of  a  tonpilz  trans¬ 
ducer  that  is  subjected  to  an  arbitrary  radiation  loading.  Since  the 
mutual  interactions  will  manifest  themselves  as  changes  in  radiation 
loading,  the  positional  dependence  can  be  accounted  for. 

The  success  of  Permutation  Search  in  solving  similar  problems 
known  as  Travelling  Salesman  Problems  indicates  it  may  be  an  effec¬ 
tive  tool  for  determining  array  positions  for  transducer  elements.  It 
is  a  straighforward  and  effective  search  technique.  For  the  problem 
of  element  placement,  it  rearranges  the  elements  in  the  array,  search¬ 
ing  for  those  permutations  which  improve  upon  the  response  of  the  cur¬ 


rent  array 


1.2  GOALS  OF  THIS  STUDY 


The  goal  of  this  thesis  is  to  develop  an  algorithm  which  deter¬ 
mines  the  placement  of  elements  in  an  array  so  that  the  impact  of  the 
transducer  element  errors  on  the  array  response  is  minimized.  This 
technique  must  be  applicable  to  all  four-fold  symmetric  arrays,  odd  and 
even.  To  accomplish  this  goal,  several  objectives  must  be  met: 

1.  Analyze  the  effects  of  tolerances  in  the  lumped-parameter  cir¬ 
cuit  values  on  an  element's  performance  characteristics. 

2.  Study  the  effects  of  mutual  interaction  between  elements  and 
the  effect  of  incidence  angle  on  the  mutual  interaction. 

3.  Adapt  the  Permutation  Search  algorithm  to  the  problem  of  ele¬ 
ment  placement,  and  determine  its  effectiveness. 

4.  Enhance  the  Kendig  Scatter  Diagram  Method  by  using  a  three  di¬ 
mensional  model  of  array  response. 

5.  Determine  the  effectiveness  of  using  the  Neo-Kendig  selection 
scheme  as  a  starting  point  for  the  Permutation  Search  algor¬ 
ithm. 

The  lumped-parameter  equivalent  circuit  and  a  study  of  the  effects 
of  parameter  tolerances  is  presented  in  Chapter  Two.  In  Chapter  Three 
is  a  discussion  of  the  effect  of  the  acoustic  coupling  between  elements 
on  an  element's  performance.  Chapter  Four  is  a  presentation  of  Permu¬ 
tation  Search  and  a  discussion  of  its  results.  The  derivation  of 


CHAPTER  TWO 


THE  TRANSDUCER  ELEMENTS 


2.1  THE  EQUIVALENT  CIRCUIT* 


In  order  to  study  the  response  characteristics  of  a  tonpilz  trans¬ 
ducer,  (shown  in  Figure  2.1),  an  equivalent  electrical  circuit  is  often 
used.  The  two  types  of  equivalent  circuits  used  are  distributed- 
parameter  and  lumped-parameter.  The  distributed-parameter  is  used  pri¬ 
marily  during  the  design  of  a  transducer.  This  study  uses  the  lumped- 
parameter  model  because  the  circuit  parameters  can  be  determined  from 
the  measurement  of  transducer  elements.  In  this  manner,  the  element 
tolerances  can  be  accounted  for  as  variations  in  the  circuit  parame¬ 
ters. 

With  the  lumped-parameter  model,  the  performance  of  a  tonpilz 
transducer  operated  in  the  vicinity  of  its  fundamental  resonance  can  be 
readily  predicted.  In  this  circuit,  the  mechanical  properties  of  the 
transducer,  mass  and  compliance,  are  represented  by  linear  electrical 
parameters,  inductance  and  capacitance.  Mechanical  resistance  is  di¬ 
rectly  analogous  to  electrical  resistance.  This  circuit  is  valid  only 
in  the  vicinity  of  the  fundamental  resonance  of  the  transducer. 

The  input/output  relationships  for  a  transducer  element  are  deter¬ 
mined  by  an  analysis  of  the  equivalent  circuit.  Shown  in  Figure  2.2 
is  the  lumped-parameter  equivalent  circuit.  For  the  receive  case,  the 
output  voltage  sensitivity,  Me ,  for  a  transducer  is  given  by 


*  The  material  in  this  section  is  a  synopsis  of  Reference  [2]. 
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Figure  2.1  A  Tonpilz  Type  Transducer 
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where  Eout  is  the  open  circuit  output  voltage  for  an  incident  sound  pres¬ 
sure  Pj>n.  From  the  equivalent  circuit,  it  can  be  shown  that 
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Equations  2-2  and  2-3,  the  input/output  relationships,  are  used  to 
determine  the  amplitude  and  phase  responses  of  each  element.  It  is  seen, 
in  Equation  2-1,  that  the  magnitude  of  the  open  circuit  voltage  is  pro¬ 
portional  to  both  the  receive  voltage  sensitivity  and  the  incident  sound 
pressure.  In  other  words,  the  output  voltage  Is  proportional  to  the  dis¬ 
placement  of  the  head  of  the  transducer,  i.e.,  amplitude. 

Included  in  the  equivalent  circuit  is  the  analytically  derived  com¬ 
plex  acoustic  radiation  impedance,  ZA,  (=  RA  +  iXA) ,  where  the  real  part, 
RA,  is  the  resistance,  and  the  imaginary  part,  XA,  the  reactance.  The 
mechanical  impedance  experienced  by  the  transducer,  Z^,  is  related  to  the 
acoustic  impedance  by  a  factor  of  1/S2:  ZA  =  Zj^/S2,  where  S  is  the  area 
of  the  face  of  the  transducer.  In  the  equivalent  circuit,  the  acoustic 
impedance  in  the  acoustic  domain  is  transformed  into  mechanical  impedance 


in  the  mechanical  domain  by  an  ideal  transformer  with  a  turns  ratio  of 
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The  input/output  relationships  show  that  a  transducer  element  will 
have  different  responses  for  different  loadings.  This  complicates  the 
array  design  process  because  an  element  placed  in  one  position  in  the 
array  will  not  experience  the  same  radiation  loading  when  placed  in  a 
different,  nonsymmetric  position.  This  is  a  result  of  the  acoustic  in¬ 
teraction  among  transducer  elements  in  an  array.  Each  element  position 
in  the  array  has  a  set  of  spatial  relationships  with  all  the  other  ele¬ 
ment  positions  in  the  array.  Those  positions  which  have  the  same  rela¬ 
tionships  with  all  the  other  array  positions  are  symmetric  positions. 

2.1.1  CIRCUIT  PARAMETERS 

A  reduced  form  of  the  lumped-parameter  circuit,  shown  in  Figure 
2.3  is  used  to  determine  the  circuit  parameters  for  each  element.  The 
relationships  between  the  two  circuits  are  shown  below. 

C0  =  CQ  C  =  CM$2  L  =  Mm/$2  R  =  Rm/®2  .  (2-4) 

In  the  reduced  form,  dielectric  losses,  R0,  are  assumed  negligible. 
Since  the  measurements  are  made  in  air,  the  radiation  impedance  is  also 
assumed  negligible.  By  measuring  Fy,  the  motional  resonance  frequency; 

and  F2,  the  half-power  frequencies;  Ct,  the  total  capacitance;  and 
Gy,  the  conductance  at  F  =  Fy,  the  circuit  parameters  of  the  reduced 
model  are  given  by 

R  =  1/Gy  C  =  Gy  (u)2  —  u)j)/wy2 

L  =  l/(C*Wy2)  c0  =  Ct  -  C,  (2-5) 

where  a  =  2*F.  The  transformation  ratio,  represents  the  coupling 
between  the  mechanical  and  electrical  domains  of  the  equivalent  circuit. 


Since  it  has  the  units  of  current /velocity  it  can  be  determined  by  the 


simultaneous  measurement  of  both  these  quantities.  A  more  detailed 
discussion  of  this  procedure  and  the  procedures  for  determining  the 
other  parameters  is  discussed  in  Reference  [2].  The  six  measured  para¬ 
meters,  Fy,  F  ,  F  ,  Gy,  Ct,  and  $ ,  are  all  that  are  needed  to  predict  a 
transducer  element's  performance  for  an  arbitrary  radiation  loading. 


2.1.2  A  MEASURED  SET 


In  order  to  simulate  realistic  sets  of  elements  to  test  the  se¬ 
lection  process,  the  relationships,  if  any,  among  the  measured  circuit 
parameters  had  to  be  determined.  A  set  of  66  tonpilz  transducers  were 
constructed  and  measured.  Correlations  and  distributions  in  measured 
data  were  tested  by  plotting  each  parameter  versus  every  other  parame¬ 
ter.  In  all  the  following  plots,  the  parameter  values  have  been  nor¬ 
malized  to  the  mean  values  of  the  measured  set. 

There  is  a  linear  correlation  between  Fy  and  both  F^  and  F2,  shown 
in  Figures  2.4  and  2.5,  respectively.  There  is  also  a  linear  correla¬ 
tion  between  Fj  and  F2,  as  seen  in  Figure  2.6.  However,  there  is  no 
correlation  between  Fy,  and  F2  -  Fj,  as  seen  in  Figure  2.7.  Neither  is 
there  a  correlation  between  Fj  or  F2  and  F2  -  Fj.  It  was  determined 


f2  =  Fy  +  (F2  "  Fj)/2 


F1  =  Fy  -  (F2  -  FL )/ 2  (2-6) 

is  a  good  representation  of  the  relationships  between  Fy,  Fj  and  F2, 
(within  .01%). 
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A  correlation  was  found  between  Gy  and  F2  -  F^,  (Figure  2.8). 


Equation  2-5,  rewritten  below,  (Equation  2-7),  shows  that  F2  -  Fj  is 


inversely  proportional  to  Gy, 


u>2  -  uj  =  C  •  ojy2/Gv 


(2-7) 


The  two  curves  shown  in  Figure  2.8  represent  the  least  squares  fit  to 


the  data  and  a  linear  approximation  to  this  curve.  The  equation  for 


the  least  sqaures  fit  is. 


F2  ~  =  .5415/Gy  +  .4402 


(2-8) 


and  the  linear  approximation  is. 


F2  -  F!  =  -.627  •  Gy  +  1.627. 


(2-9) 


This  relationship  is  shown  in  Figure  2.9,  a  plot  of  conductance,  G, 


versus  frequency.  As  the  value  of  G  at  resonance,  Gy,  increases,  the 
reson’  peak  becomes  narrower.  Therefore,  an  increase  in  Gy  produces 


a  de-  se  in  F2  -  F^. 


These  plots,  along  with  Figures  2.10  and  2.11,  have  shown  that  Gy, 
Fy,  and  Ct  have  approximately  Gaussian  distributions.  Additionally,  it 


is  seen  that  $  has  a  random  distribution.  This  is  not  surprising  since 


$  represents  the  coupling  between  the  electrical  and  mechanical  domains 


of  the  transducer.  Random  construction  errors  show  up  as  random  fluc¬ 


tuations  in  # . 


From  the  relationships  expressed  in  Equations  2-6  and  2-7,  it  is 


evident  that  instead  of  the  six  measured  parameters,  there  are  only 


four  independent  parameters  required  to  describe  the  performance  of  a 


transducer:  Gy,  Fy,  $,  and  Ct. 
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A  simulated  set  of  elements  was  generated  using  the  results  of 
this  study.  Fy,  Gy,  and  Ct  were  generated  by  a  Gaussian  number  gener¬ 
ator  with  standard  deviations,  a,  expressed  as  percentages  of  the  mea¬ 
sured  means,  of  0.9%,  17.4%,  and  3.9%,  respectively.  A  random  number 
generator  and  a  ±12%  distribution  was  used  to  generate  values  for 
The  resulting  selection  pool,  when  represented  in  a  scatter  plot  of  am¬ 
plitude  vs.  phase,  was  indistinguishable  from  the  real  selection  pool. 
The  actual  range  of  values  for  each  parameter  is  ±4%  for  Fy,  ±45%  for 
Gy,  ±11%  for  Ct,  and  ±13%  for  $.  These  ranges,  in  turn,  produce  ampli¬ 
tude  and  phase  tolerances  of  ±9%  and  ±11°,  respectively. 

2.2  A  TOLERANCE  STUDY 

One  of  the  obvious  results  of  this  thesis  is  that  better  response 
patterns  can  be  obtained  if  the  amplitude  and  phase  responses  of  the 
elements  do  not  vary  greatly  from  one  another.  Thus,  it  is  advanta¬ 
geous  to  remove  from  the  selection  pool  those  elements  whose  amplitude 
and  phase  responses  do  not  fall  within  a  specified  tolerance  limit  of 
the  mean  values. 

Calculating  the  amplitude  and  phase  responses  of  each  element  in 
every  array  position  is  a  lengthy  process.  If  relationships  between 
circuit  parameter  values  and  element  performance  characteristics  can  be 
found,  these  relationships  could  be  used  as  a  criterion  for  element 
suitability,  thereby  avoiding  lengthy  calculations. 


Using  the  input/output  relationships  of  the  equivalent  circuit,  it 
is  possible  to  estimate  the  effects  of  extreme  values  for  each  para¬ 
meter,  and  establish  acceptable  parameter  value  limits.  The  amplitude 
and  phase  responses  of  a  nonspecific  transducer  are  calculated  as  one 
parameter  value  is  changed  a  specified  percentage  above  and  below  a 
measured  mean  value  and  the  other  parameters  are  held  constant.  The 
mean  value  is  obtained  from  the  measurement  of  a  set  of  elements.  The 
specified  percentage  is  based  on  the  percent  differences  between  the 
mean  value  and  the  maximum  and  minimum  values  of  the  measured  set.  To 
see  how  the  parameter  values  interact,  three  lines  are  plotted  on  each 
graph.  Each  line  represents  different  fixed  values  of  an  additional 
parameter.  This  discussion  is  limited  to  the  effects  on  an  element's 
open  circuit  voltage  amplitude  response. 

Those  parameters  for  which  the  amplitude,  (and  phase),  responses 
of  the  elements  do  not  change  significantly  over  their  entire  ranges 
are  $  and  Gy,  (Figure  2.12).  The  amplitude  response  of  an  element  does 
not  vary  significantly  for  a  ±30%  variation  in  $,  or  a  ±80%  in  Gy; 
large  changes  in  Gy  or  4  will  not  produce  large  changes  in  amplitude. 
In  contrast,  Figure  2.13  shows  the  effect  of  a  ±30%  variation  in  Ct, 
and  ±20%  variation  in  Fy.  Large  changes  in  either  of  these  parameters 
will  produce  a  large  change  in  amplitude  response. 

Although  the  four  circuit  parameters  are  independent,  one  para¬ 
meter  can  influence  the  effect  another  will  have  on  an  element's  re¬ 
sponse.  Figure  2.14  shows  how  the  amplitude  response  of  an  element 
changes  as  a  function  of  Gy,  with  Fy  as  the  additional  parameter.  The 
relationships  mentioned  in  the  previous  paragraph  are  clearly  visible. 


Note  that  near  the  lower  extreme  of  Gy  there  is  a  sudden  change  in  slope 
for  the  curve  indicating  Fy  at  five  percent  below  its  mean  value.  This 
indicates  that  the  effect  of  Gy  on  amplitude  is  somewhat  coupled  to  the 
value  of  Fy 

Figure  2.15  is  a  similar  plot,  but  Fy  is  allowed  to  vary  over  a 
larger  range.  As  seen,  for  Fy  at  twenty  percent  below  its  measured  mean 
value,  the  coupling  between  Gy  and  Fy  is  more  pronounced.  This  results 
in  a  reduced  range  of  acceptable  values  for  Gy.  Therefore,  a  greater 
tolerance  in  Fy  reduces  the  range  of  acceptable  values  of  Gy.  Trans¬ 
ducers  with  Fy  and  Gy  values  in  this  range  may  be  unsuitable  for  the  ar¬ 
ray. 

Although  a  ±20%  tolerance  on  Fy  violates  the  lumped-parameter  equi¬ 
valent  circuit  assumption,  this  plot  is  a  valid  representation  of  the 
coupling  of  Fy  and  Gy  And,  even  though  the  limits  on  the  parameter  val¬ 
ues  in  these  plots  far  exceeds  those  found  in  the  measured  set,  these 
plots  do  show  the  relative  importance  of  each  parameter.  The  results  in¬ 
dicate  that  the  values  of  ♦  and  Gy  are  not  as  influential  on  the  response 
characteristics  as  are  Fy  and  Ct.  It  follows  that  greater  uniformity  of 
element  performance  characteristics  can  be  achieved  by  closely  control¬ 
ling  the  tolerances  on  Fy  and  Ct. 

Since  only  one  set  of  elements  was  measured,  it  is  not  possible  to 
set  absolute  or  relative  limits  on  the  parameter  values  for  other  sets 
based  on  this  set.  However,  this  method  has  shown  the  effects  each  of 
the  circuit  parameters  has  on  an  element's  amplitude  and  phase  response. 


and  on  each  other 


2.3  ARRAY  RESPONSE  AND  ELEMENT  TOLERANCES 


Prior  to  this  study,  a  rule  of  thumb  relationship  between  array 
response  and  constituent  element  tolerances  was  to  set  the  tolerances 
on  amplitude  and  phase  equal  to  the  desired  side  lobe  level.  For  exam¬ 
ple,  to  achieve  a  -40  dB  side  lobe  level,  a  .01  tolerance  in  amplitude 
and  phase  is  required.  This  translates  into  a  - 1 %  tolerance  in  ampli¬ 
tude  and  a  ±0.6°  tolerance  in  phase.  These  tolerances  are  very  strin¬ 
gent  ,'.rid  very  difficult  to  obtain. 

The  selection  process  presented  in  this  thesis  changes  the  above 
mentioned  tolerances  by  an  order  of  magnitude.  Elements  with  ±9%  am¬ 
plitude  and  ±11°  phase  tolerances,  when  randomly  placed  in  an  array 
produce  -30  dB  side  lobe  levels,  as  predicted  by  the  rule  of  thumb  re¬ 
lationship.  However,  using  the  same  elements,  the  selection  process 
can  produce  arrays  with  side  lobe  levels  that  are  within  1  dB  of  the 
design  level  of  -40  dB  side  lobes,  (see  Chapter  Six). 


CHAPTER  THREE 


IMPEDANCE  AND  INTERACTION 


3.1  INTRODUCTION 

Mechanical  impedance,  Zf^,  is  defined  as  the  ratio  of  the  force  ap¬ 
plied  to  an  object  and  the  resulting  velocity  at  the  point  of  applica¬ 
tion.  In  general,  impedance  is  a  complex  quantity.  The  real  part,  re¬ 
ferred  to  as  resistance,  is  related  to  the  power  lost  in  the  system. 
The  imaginary  part,  reactance,  represents  energy  which  is  stored  in  the 
system. 

In  acoustics,  a  more  useful  quantity  is  the  acoustic  impedance, 
ZA*  It  is  defined  as  the  ratio  of  the  pressure  on  a  surface  and  the 
volume  velocity  of  the  medium  which  results  from  the  motion  of  the  sur¬ 
face.  Since  pressure  and  force  are  related  by  a  factor  of  area,  as  are 
volume  velocity  and  particle  velocity,  mechanical  impedance  and  acous¬ 
tic  impedance  are  related,  Z^  =  •  Z A,  where  S  is  the  area  of  the 
face  of  the  transducer. 

For  a  plane  wave  incident  on  a  piston,  the  force  exerted  on  the 
piston  is  obtained  by  averaging  the  pressure  over  the  surface  area  of 
the  piston.  In  this  way,  the  acoustic  impedance  is  transformed  into 
mechanical  impedance.  In  the  equivalent  circuit,  the  coupling  of  the 
mechanical  domain  with  the  acoustic  domain  is  represented  by  an  ideal 
transformer  with  a  turns  ratio  equal  to  the  area  of  the  transducer 


face. 
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3.1.1  SELF  IMPEDANCE 

The  head  of  a  tonpilz  transducer  can  be  modelled  as  a  rigid  square 
piston  in  an  infinite  baffle.  For  a  piston  which  is  small  compared  to 
wavelength,  the  pressure  on  the  face  of  the  piston  will  vary  with  posi¬ 
tion.  This  being  the  case,  the  acoustic  impedance  of  the  transducer  is 
no  longer  a  simple  ratio. 

Using  the  Green's  function  for  radiation  from  a  planar  surface, 
Morse  and  Ingard  [3]  derived  an  expression  for  the  acoustic  radiation 
impedance  of  a  square  piston  of  side  a.  The  expression  is  given  below. 


Z(w)  =  pea  [6  (ka)  +  jx  (ka)] 
o  o 


(3-1) 


m 


where 


.  *12  „ 

X  =  M  (ka)  =  —  /  sin((ka)cos(q))sin  (q)dq 
°  0 


0  =  1  -  2J,(ka)/ka 

o  I 


•y&l 


3.1.2  MUTUAL  IMPEDANCE 

For  the  two  square  pistons  shown  in  Figure  3.1,  the  pressure  on 
the  face  of  one  piston  is  a  function  of  its  own  interaction  with  the 
medium  and  of  the  pressure  field  of  the  other  transducer.  Simultane¬ 
ously,  the  other  piston  is  experiencing  the  same  effects  as  a  result  of 
the  pressure  from  the  first  piston.  The  motion  of  each  piston  is, 
therefore,  coupled  with  the  other. 
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As  expressed  below,  the  total  impedance  on  one  transducer  is  the 
combination  of  its  self  impedance, ,  and  the  mutual  impedance  with 
the  other  transducer. 


ZT  ’  Zs  +  <W  •  *1  2 


(3-2) 


1  2 

Z  is  called  the  mutual  impedance  coefficient.  v„/v,  represents  the 
m  2  1 

of  the  velocity  amplitudes  of  the  two  transducers.  The  superscripts 
refer  to  the  elements  involved. 

For  an  array  of  M  transducers,  the  total  loading  on  one  element  is 
a  function  of  its  self  impedance  and  the  weighted  sum  of  the  mutual  im¬ 
pedances  with  all  the  other  elements. 


2 

T 


M 


+  Z  (v  /v  .  ) 
i  =  l  J 


(3-3) 


If  all  the  can  be  determined,  as  well  as  Z ^  ,  it  would  be  possible 
m  s 

to  determine  the  acoustic  radiation  loading  on  every  element  in  an  ar¬ 
ray.  The  ratio,  Vj/vj,  representing  the  ratio  of  velocity  amplitudes, 
is  equivalently  a  ratio  of  shading  coefficients. 

Using  the  expression  for  self  impedance,  Equation  3-1,  and  an  ex¬ 
pression  for  the  mutual  impedance.  Reference  [4],  Schafer  [2]  has  de¬ 
veloped  a  computer  routine,  "MUTUAL",  to  compute  the  radiation  loadings 
on  all  the  positions  of  an  8  by  8  array.  For  this  study,  this  program 


has  been  expanded  to  handle  all  types  of  four-fold  symmetric  arrays  up 
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3.2  INTERACTION  FOR  THE  TRANSMIT  AND  RECEIVE  CASES 


Since  the  mutual  interaction  depends  on  a  ratio  of  the  velocities 
of  two  elements  vibrating  in  phase,  a  change  in  shading  values  will  af¬ 
fect  this  interaction.  When  an  array  is  used  in  the  transmit  mode,  the 
elements'  amplitudes  are  shaded  in  a  manner  to  produce  a  specified  di¬ 
rectional  response.  Each  set  of  shading  values  corresponding  to  a  dif¬ 
ferent  beam  pattern  will  subject  the  array  elements  to  different  mutual 
interaction  effects. 

The  receive  case  is  different.  Although  shading  values  are  still 
used  in  beamforming,  they  have  no  effect  on  the  acoustic  interaction 
between  elements.  The  shading  values  are  applied  to  the  output  elec¬ 
trical  signals  and  therefore  do  not  influence  the  impedance  on  or  the 
motion  of  the  transducers. 


3.2.1  THE  INFLUENCE  OF  THE  INCIDENT  ANGLE 


Unfortunately,  the  receive  case  can  not  be  simply  modelled  by 
equating  it  to  the  transmit  case  with  unity  shading.  For  a  plane  wave 
with  an  incident  angle  9+0°,  the  elements  separated  by  d  =  A/2  will 
not  vibrate  in  phase.  If  the  elements  are  modelled  as  point  sources 
placed  at  their  respective  centers,  the  responses  from  two  elements 
separated  by  a  distance  d  will  be  out  of  phase  by  k*d*sin9,  (see  Figure 
3.2). 

If  the  phase  relationship  between  two  elements  changes,  so  will 
their  interaction.  The  interaction  is  reflected  in  the  radiation  re¬ 
sistance  and  the  power  radiated.  To  see  this,  consider  two  sources 
separated  by  a  distance  d.  If  these  two  elements  vibrate  180°  out  of 
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phase,  as  in  a  dipole,  the  total  power  radiated  is,  for  k'd  =  1,  one 
third  of  that  for  a  single  source.  As  k*d  changes,  so  does  the  amount 
of  power  radiated.  When  two  elements  vibrate  90°  out  of  phase,  the  to¬ 
tal  power  radiated  is  twice  that  of  a  single  source.  Similarly,  for 
two  sources  vibrating  in  phase,  the  total  power  radiated  is  four  times 
that  of  a  single  source  [5]. 

It  is  easy  to  extend  this  example  to  the  receive  case.  The  result 
is  that  as  the  incident  angle  changes,  so  do  the  phase  relationships, 
and  the  elements'  mutual  interactions.  Consequently,  each  element's 
performance  characteristics  will  change  with  the  incident  angle. 

This  is  disastrous  news  for  the  problem  of  placing  elements  in  an 
array.  If  the  phase  relationships  are  determined  for  a  particular  in¬ 
cident  angle,  the  radiation  loadings  and  element  responses  can  be  de¬ 
termined.  An  array  can  then  be  designed  which  has  minimal  error  for 
that  incident  angle.  However,  the  errors  at  other  incident 
angles  are  not  accounted  for  by  this  element  placement. 


3.3  TRADITIONAL  AND  NEW  ASSUMPTIONS 

Before  the  inclusion  of  mutual  interaction  effects  became  fea¬ 
sible,  it  was  assumed  that  the  entire  array  experienced  Pc  loading. 
Since  the  results  based  on  this  assumption  have  been  acceptable,  the 
mutual  interaction  effects  were  considered  second  order  effects. 
Through  the  use  of  "MUTUAL",  it  is  shown  that  mutual  interaction  ef¬ 
fects  on  individual  elements  are  not  second  order. 
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In  the  transmit  case,  the  loading  on  the  centermost  elements  dif¬ 
fers  from  the  loading  on  a  single  square  element,  as  computed  using 
Equation  3-1,  by  6%  in  the  resistive  component  and  -83%,  reactive. 
This  error  increases  to  18%,  resistive,  and  -97%,  reactive  for  elements 
farthest  from  the  center.  In  the  receive  case,  for  a  0°  incident  an¬ 
gle,  the  errors  are  9%,  resistive  and  -85%,  reactive  for  centermost  and 
4.5%,  resistive  and  -52%,  reactive  for  outermost. 

These  errors  are  far  from  what  may  reasonably  be  considered  second 
order.  Although  the  loading  on  every  element  is  dramatically  different 
from  pc  loading,  the  loading  on  the  entire  array,  with  interaction  ef¬ 
fects  included,  can  still  be  reasonably  approximated  by  pc  loading, 
[6].  The  larger  the  array,  the  more  accurate  is  the  approximation. 

For  the  purpose  of  this  study,  it  is  assumed  that  for  plane  waves 
incident  at  angles  *  0°,  the  radiation  loadings  on  the  elements  of  an 
array  do  not  differ  significantly  from  the  loadings  for  0°  incidence. 
Thus,  the  radiation  loadings  for  all  the  incident  angles  are  assumed  to 
be  the  same  as  the  loadings  for  0°  incidence. 

The  following  section  illustrates  how  changes  in  radiation  loading 
do  not  significantly  affect  element  performance.  Therefore,  the  small 
changes  in  radiation  loading  which  occur  over  different  incident  angles 
will  not  produce  significant  changes  in  element  response. 


3.4  THE  EFFECT  OF  ELEMENT  TOLERANCES  ON  MUTUAL  INTERACTION 


In  order  to  determine  the  radiation  loading  on  every  position  in 
an  array,  it  is  necessary  to  know  the  shading  values  for  each  position. 
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These  radiation  loadings  are  then  used  to  determine  the  elements'  per¬ 
formance  tolerances.  However,  the  use  of  ideal  shading  values  assumes 
that  each  shaded  element  is  identical.  But  this  is  not  true!  It  is 
important  to  ask,  and  answer,  what  effects  do  element  tolerances  have 
on  the  radiation  loading? 

An  array  was  designed  using  the  ideal  radiation  loadings  to  cal¬ 
culate  element  tolerances.  The  normalized  mean  amplitude  values  of  the 
elements  in  each  symmetric  array  position  were  then  multiplied  by  the 
shading  values  for  each  position  to  produce  new  shading  values  for  each 
position.  It  was  assumed  that  the  phase  errors  of  each  element  are 
small  and  were  not  included.  These  new  shading  values  were  then  used 
to  re-evaluate  the  radiation  loadings  for  each  array  position. 

For  those  elements  which  comprise  this  array,  new  amplitude  and 
phase  characteristics  were  determined.  A  new  response  pattern  for  this 
array  was  then  determined  using  the  newly  computed  element  characteris¬ 
tics.  The  resulting  response  pattern  did  not  exhibit  significant 
changes  from  the  original  pattern.  The  new  radiation  loading  that  the 
elements  experienced  due  to  their  neighbors'  new  tolerances  was  also 
computed. 

The  first  recalculation  of  the  radiation  loadings  on  each  position 
showed  changes  ranging  from  0%  to  3%  in  resistance  and  1.5%  to  90%  in 
reactance.  Yet  these  changes  resulted  in  a  maximum  change  of  only  3% 
in  amplitude  and  0.5°  in  phase.  It  should  be  noted  that  the  greatest 
changes  occur  for  radiation  loadings  corresponding  to  the  array  posi¬ 
tions  farthest  from  the  array  center.  The  errors  for  the  closer 
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positions  are  an  order  of  magnitude  smaller 


The  second  recalculation 


of  radiation  loadings  produced  even  smaller  changes.  This  translates 
into  smaller  errors  in  the  amplitude  and  phase. 

Additionally,  these  new  performance  characteristics  were  used  to 
redesign  the  array  from  the  original  selection  pool.  It  was  assumed 
that  the  errors  in  each  symmetric  array  position  would  not  change  sig¬ 
nificantly  from  those  of  the  first  array.  (This  was  tested  and  shown 
to  be  valid.)  A  comparison  of  this  new  response  pattern  with  the  first 
also  showed  no  significant  degradation  or  improvement.  The  normalized 
mean  amplitude  values  for  each  position  were  used  to  recalculate  the 
amplitude  and  phase  characteristics  of  each  element.  These  new  charac¬ 
teristics  were  not  significantly  different  from  the  original  values. 

In  conclusion,  the  errors  introduced  into  the  radiation  loadings 
by  element  tolerances  are  small  enough  to  be  neglected. 


CHAPTER  FOUR 


PERMUTATION  SEARCH 


4.1  INTRODUCTION 

The  process  of  placing  in  an  array  elements  which  have  different 
tolerances  involves  combining  the  elements  in  such  a  way  as  to  minimize 
the  total  impact  of  their  errors  on  the  directional  beam  pattern  of  the 
array.  This  process  resembles  a  class  of  problems  in  industrial  engi¬ 
neering  known  as  the  Travelling  Salesman  Problem.  Travelling  Salesman 
Problems  involve  ordering  a  set  of  elements  in  different  ways  in  an  at¬ 
tempt  to  meet  a  certain  objective.  The  elements  could  represent  the 
stops  a  salesman  has  to  make  during  a  trip  and  the  objective  could  be 
to  minimize  the  number  of  miles  travelled.  Or,  in  the  case  of  array 
design,  each  element  would  represent  a  transducer  element  placed  in  an 
array  position,  and  the  objective  would  be  to  minimize  the  error  impact 
on  the  array  response. 

For  small  problems,  i.e.,  those  involving  a  small  number  of  ele¬ 
ments,  there  are  solution  techniques  which  will  give  a  globally  optimum 
solution.  But  with  large  problems,  the  amount  of  computer  time  these 
methods  require  becomes  unreasonable.  For  an  array  of  36  positions  and 
36  elements,  there  are  36!  possible  permutations  of  the  elements. 
There  are  other  methods  of  solution  which  may  not  produce  a  globally 
optimum  solution,  but  do  offer  a  considerable  reduction  in  necessary 
computer  time.  One  such  method  is  called  Permutation  Search. 


4.1.1  PROBLEM  DEFINITION 


A  permutation,  denoted  [Pn] ,  is  any  ordering  of  a  set  of  elements. 
For  each  permutation,  there  is  an  associated  objective  function,  E(Pn). 
Additionally,  there  may  be  constraints  on  how  these  elements  may  be  ar¬ 
ranged  in  the  permutation,  denoted  as  the  set  of  feasible  permutations, 
G.  The  problem  solved  by  Permutation  Search  can  be  formulated  as 


Find  [Pn]  to  minimize 
E(Pn),  such  that 
Pn  EG. 


4.1.2  THE  ELEMENT  PLACEMENT  PROBLEM 


(4-1) 


For  the  process  of  creating  an  array  with  M  positions,  using  a  set 
of  N  transducer  elements,  where  N  >  M,  the  problem  formulation  is  as 


follows. 


Let  X(m,n)  be  the  decision  variable  associated  with  element  m  and 
position  n,  such  that 


X(m,n)  =1  if  element  m  is  in  position  n 
=  0  otherwise. 


(4-2) 


There  are  two  constraints  on  the  placement  of  the  elements.  The  first 
is  that  every  array  position  must  be  filled  and  by  one  element  only. 
For  those  humans  who  know  that  two  elements  may  not  occupy  the  same  po¬ 
sition  simultaneously,  the  second  part  of  this  constraint  may  seem  un¬ 
necessary.  However,  computers  have  no  way  of  knowing  this  unless  they 
are  told  so. 
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N 

£  X(m,n)  =  1 
m=l 


n  =  1  , 2  , . . .  ,M 


(4-3) 

The  second  constraint  is  that  each  element  may  either  be  assigned  to 
one  array  position,  or  unassigned. 


M 

£  X(m,n)  <  1 
n=l 


(4-4) 


The  array  response,  R(  0,4)),  at  incident  angles  9  and  <t>,  is 


R(0,4>)  =  20  log 


where 


M  N 

[  £  £  X(m,n)  A(m,n) • exp (i [kun  +  a(m,n)]}J 

m=l  n=l 

(4-5) 


X(m,n)  =  Decision  variable 
A(m,n)  =  S(n)  •  Real (Elemnt (m,n) } 

S(n)  =  Shading  value  at  position  n 

N  =  Number  of  array  positions 

M  =  Number  of  elements  in  the  selection  pool 

a(m,n)  =  Irnaginary{Elemnt(m,n) 1 

un  =  k[C(n)*sin0  cos<t>  +  D(n)*sin9  sim}>] 

C'n)  =  x  coordinate  of  position  n 

D(n)  =  y  coordinate  of  position  n 

Elemnt(m,n)  =  Amplitude,  (real)  and  Phase,  (imaginary) 
response  of  element  m  at  position  n. 


The  objective  is  to  minimize  the  total  error,  defined  as  the  sum  over 
all  constrained  angles  of  the  error  at  each  angle: 


E  (X(m,n)  }  =  I  (L(  9 ,  <f> )  -  R(0,j>)}  •  e, 
B 


(4-6) 


where 


L(0,4>)  =  Desired  level  at  angles  9  and 

B  =  The  set  of  all  constrained  angles 
e  =  1  if  L(9,<*>)  >  R(  9, 6) 

=  0  if  L(  9 ,  <f>)  <  R(9,4>). 
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4.2  THE  SEARCH  PROCESS* 


Permutation  Search  is  a  process  which  searches  for  a  locally  opti¬ 
mum  permutation.  It  looks  in  the  neighborhood  of  an  initial  base  per¬ 
mutation  for  permutations  which  improve  the  objective  function.  The 
neighborhood  of  a  permutation  is  defined  as  the  set  of  permutations 
which  can  be  obtained  through  a  series  of  element  exchanges.  A  base 
permutation  is  any  permutation  from  which  a  new  permutation  is  formed. 
The  initial  base  permutation,  or  the  first  permutation  in  a  series  of 
permutations,  is  denoted  [P0].  For  example,  a  permutation,  [P*-],  of  5 
e lements , 

[P1 ]  =43251 

is  in  the  neighborhood  of  the  permutation 
[P°]  =54321 

since  it  can  be  reached  from  [p°]  through  a  sequence  of  one  for  one  ad¬ 
jacent  element  exchanges,  viz., 

[P°]  =54321 
[  P1  ]  =  4  5  3  2  1 
[P2]  =  4  3  5  2  1 
I P3 1  =43251. 

When  the  final  element  in  the  permutation  is  exchanged,  the  next 
exchange  involves  the  first  two  elements  of  the  existing  permutation. 
In  thr  ample  above,  the  next  two  permutations  would  be 
fP4]  =432  1  5 
[P5]  =34215. 

*  The  Permutation  Search  orocedure  presented  in  this  section  is  a 
synopsis  of  reference  [7] 
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Exchanges  can  be  the  interchanging  of  two  elements,  or  shifts  of 
large  blocks  of  elements.  The  nature  of  the  problem  dictates  the  ap¬ 
propriate  type  of  exchange.  For  the  transducer  array,  one  for  one  ex¬ 
changes  of  adjacent  elements  are  performed. 

4.2.1  EVALUATING  THE  PERMUTATIONS 

Since  Permutation  Search  looks  for  those  permutations  which  im¬ 
prove  the  value  of  the  objective  function,  E(Pn),  E  must  be  determined 
for  every  permutation.  It  is  obvious  that  not  every  exchange  will  re¬ 
sult  in  a  permutation  with  an  improved  objective  function.  If  an  im¬ 
provement  is  found,  then  the  next  exchange  is  performed  on  the  newly 
achieved  permutation,  at  the  position  adjacent  to  the  previous  ex¬ 
change.  If  there  is  no  improvement,  the  next  exchange  is  performed  on 
the  previous  permutation,  also  at  the  adjacent  position.  Any  per¬ 
mutation  obtained  as  a  result  of  an  element  exchange  performed  on  a 
base  permutation  is  called  a  trial  permutation.  If  this  trial  per¬ 
mutation  improves  the  value  of  the  objective  function  it  becomes  the 
current  base  permutation.  Therefore,  a  sequence  of  bases  represents  a 
continuous  improvement  in  the  objective  function.  This  sequence  will 
converge  to  a  locally  optimum  solution  in  the  neighborhood  of  the  ini¬ 
tial  base  permutation. 

4.3  ARitAY  DESIGN  SIMPLIFICATION 

The  Permutation  Search  formulation  of  element  placement  problems 
can  be  set  up  in  such  a  way  as  to  sidestep  the  two  constraints.  Let 
the  permutation  array,  (note,  not  permutation),  be  a  vector  with  as 


many  components,  or  slots,  as  there  are  transducer  elements  in  the  se¬ 
lection  pool.  Each  slot  will  represent  either  a  position  in  the  trans¬ 
ducer  array,  or  a  position  on  a  workbench,  i.e.,  an  unassigned  posi¬ 
tion.  A  permutation  consists  of  each  transducer  element  assigned  to 
only  one  slot  in  the  permutation  array.  Thus,  every  slot  in  the  permu¬ 
tation  array  is  filled.  This  arrangement  automatically  ensures  that 
any  permutation  in  the  neighborhood  of  the  first  permutation  satisfies 
the  constraints. 

Finding  R( © , 4> )  can  also  be  simplified  through  the  implementation 
of  the  permutation  array.  Let  there  be  a  function  assigned  to  each 
slot  of  the  permutation  array.  This  function  operates  on  the  amplitude 
and  phase  characteristics  of  the  element  assigned  to  that  slot  in  a 
manner  which  is  unique  to  the  array  position  that  the  slot  represents. 
The  uniqueness  of  this  function  is  derived  from  the  x  and  y  coordinates 
of  each  position  as  well  as  the  shading  value  at  the  position.  R  is 
found  by  summing  the  value  of  this  function  over  all  the  slots.  Thus, 
the  problem  reduces  to 

Given  [P° ] , 

Find  [Pn]  to  minimize  (4-7) 

E(Pn) . 

4.3.1  PERMSER 

The  application  of  Permutation  Search  to  the  problem  of  element 
placement  in  arrays  is  best  described  by  an  example.  This  example  dem¬ 
onstrates  how  it  is  implemented  in  the  FORTRAN  routine,  "PERMSER",  (see 
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Appendix  D  for  a  listing  and  section  4. A  for  a  discussion  of  this  pro¬ 
gram).  There  are  seven  elements  in  the  selection  pool,  numbered  1 
through  7.  There  are,  therefore,  seven  slots  in  the  permutation  array, 


[P  1  = 


(4-8) 


1  2  3  4  5 


6  7 


The  "transducer  array”  to  be  designed  has  five  positions,  represented 
by  the  first  five  slots  in  the  permutation  array.  "Workbench"  posi¬ 
tions,  slots  6  and  7,  indicate  an  element  unassigned  to  an  array  posi¬ 
tion.  Normally,  transducer  elements  are  assigned  a  number  to  distin¬ 
guish  them.  For  this  example,  the  “response”  of  each  "transducer  ele¬ 
ment"  will  be  represented  by  its  assigned  number.  The  unique  function 
at  each  slot,  representing  the  contribution  to  the  array  response  of 
each  element,  will  be  the  product  of  the  slot  index  and  the  "response” 
of  each  "transducer  element".  "Workbench"  elements  do  not  contribute 
to  the  response  of  an  array. 

The  response  of  the  array  is  the  sum  of  the  values  of  the  function 
at  each  slot  over  all  the  "transducer  array  position”  slots.  For  this 
example  the  objective  will  be  to  maximize  E(Pn). 

Starting  with  the  initial  base  permutation, 


E[P°)  =  64 


[P°]  =14265/37 


the  first  exchange  gives  a  trial  perrauation, 


[  Pt  ]“iI165  /  37_ 


with  an  objective  function  E[Pt]  =  61.  This  permutation  does  not  re¬ 
sult  in  an  improved  objective  function,  so  the  next  exchange  is  per¬ 
formed  at  the  adjacent  permutation  array  slot  in  [P°], 
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[pt!  =i24M  /  3  7_ 


The  objective  function  is  E[Pt]  =  66.  Since  this  is  an  improvement 


over  the  value  of  the  objective  of  the  current  base  point,  this  trial 


permutation  becomes  the  current  base  point, 


[P1 )  *12465/37 


EfP1]  =  66. 


Exchanging  cyclically  through  the  permutation  array,  the  process  con¬ 


verges  to  a  locally  optimum  solution.  The  sequence  of  successive  bases 


is  given  below. 


[P°]  *  1  4  2  6  5  /  3  7 


E[P°]  =  64 


[pll  =11  4  6  5  /  3  1_ 


E [ P 1 ]  =  66 


[P2]  =1  2456/  37 


E[P2]  =  67 


A  local  optimum  has  been  reached.  At  this  point,  "PERMSER" ,  will 


initiate  what  is  called  a  new  start.  A  new  permutation  is  generated  by 


randomly  rearranging  the  unassigned  elements,  and  leaving  the  elements 


which  are  assigned  to  "transducer  array  position"  slots  where  they  are. 


It  is  called  a  new  start  because  the  new  permutation  represents  a  new. 


untested  neighborhood  and  retains  the  same  level  of  achievement  of  the 


objective  function. 


The  new  start,  shown  below,  has  the  same  objective  function  value 


as  the  current  base  permutation,  but  is  a  different  permutation. 


[P3]  *12456/73 


E ( P3 ]  =  67 
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From  here,  there  is  only  one  more  base  permutation  that  will  improve 
the  objective  function.  The  new  local  optimum  is 


[P 4 1=11^5  7^/  6  3 


The  global  optimum. 


[Pn]  =34567/12 


E[P4]  =  72. 


E[Pn]  =  85 


is  not  attainable  in  either  of  the  two  neighborhoods  searched. 

In  this  example,  the  elements  in  the  permutation  array  slots  rep¬ 
resenting  the  workbench  were  never  in  a  "transducer  array  position" 
slot  for  a  base  permutation,  (except  for  the  final  permutation).  This 
particular  type  of  permutation  array  is  very  restrictive  in  the  sense 
that  there  is  only  one  entry  point  into  the  "transducer  array  position" 
slots  for  elements  that  are  in  "workbench"  slots.  A  much  more  effect¬ 
ive  permutation  array,  which  has  many  entry  points,  is  illustrated  be- 


*  * 
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(4-9) 


For  this  permutation  array,  "transducer  array  positions"  are  represent¬ 
ed  by  slots  1,  3,  4,  6,  7,  and  the  "workbench”  slots  are  2,  and  5. 
This  arrangement  gives  the  "transducer  elements"  in  the  "workbench" 
slots  a  greater  chance  of  being  exchanged  into  a  "transducer  array  po¬ 
sition"  slot. 


The  emphasis  on  the  ease  of  getting  elements  from  the  workbench 
into  the  transducer  array  is  justified.  Due  to  the  way  the  responses 
of  elements  in  symmetric  array  positions  interact,  it  is  important  that 
as  many  new  combinations  as  possible  be  tried,  (see  section  5.4.1). 
This  requires  that  elements  initially  on  the  workbench  be  included  as 
often  as  elements  initially  in  the  transducer  array. 


4.4  TEST  RESULTS 


"PERMSER",  the  FORTRAN  implementation  of  Permutation  Search,  has 
been  successfully  used  to  place  elements  in  both  3  by  3  and  6  by  6  ar¬ 
rays  from  scratch,  as  well  as  7  by  7 ,  8  by  8,  and  9  by  9  arrays  with 
non-random  starting  points.  It  can  handle  any  four-fold  symmetric  ar¬ 
ray  as  well  as  most  other  types  of  arrays.  There  are  several  param¬ 
eters  which  are  used  to  control  the  Permutation  Search  algorithm,  each 
with  different  effects  on  the  outcome.  They  can  influence  run  time, 
quality  of  result,  both  or  neither.  The  following  parameters  are  those 
which  are  not  directly  specified  by  the  problem. 


4.4.1  ZLOLIM 


The  computer  variable  ZLOLIM  represents  the  desired  or  acceptable 
amount  of  error  in  the  response  pattern  of  the  array  being  designed.  A 
comparison  of  ZLOLIM  with  the  value  of  the  objective  function  at  each 
base  permutation  determines  if  the  desired  error  has  been  achieved. 
ZLOLIM  controls  both  the  run  time  and  the  quality  of  results.  By  set¬ 
ting  ZLOLIM  very  low,  the  run  time  increases,  and  the  quality  of 
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results  will  usually  improve.  The  improvement  in  results  is  a  reflec¬ 
tion  of  the  number  and  size  of  the  neighborhoods  searched.  If  more 
time  is  used,  more  neighborhoods  can  be  searched. 

The  value  of  ZLOLIM  can  be  determined  by  using,  in  the  YVALUE  rou¬ 
tine  of  "PERMSER",  the  parameters  of  an  array  whose  response  is  already 
known.  An  evaluation  of  how  well  this  array  meets  the  design  specifi¬ 
cations  at  the  constrained  angles  will  give  a  ballpark  figure  for 
ZLOLIM. 

An  alternative  use  of  ZLOLIM  is  to  set  it  equal  to  zero.  The  run 
time  is  then  controlled  by  other  search  routine  control  parameters. 
This  method  results  in  a  loss  of  control  over  the  amount  of  acceptable 
error,  but  eliminates  the  need  to  determine  ZLOLIM  in  advance.  This  is 
the  first  of  the  two  control  parameters  which  are  capable  of  stopping 
the  execution  of  the  Permutation  Search  algorithm. 


4.4.2  CONSTRAINED  ANGLES 


Another  parameter  which  controls  the  Permutation  Search  algorithm 
is  the  number  and  location  of  the  constrained  angles.  As  the  number  of 
constrained  angles  increases,  so  does  the  run  time.  The  only  limits  on 
the  maximum  number  of  constraint  points  are  those  imposed  by  the  com¬ 
puter  facilities  and  the  amount  of  time  allowed.  However,  there  are 
limits  to  the  fewest  number  of  constraints  which  can  be  specified.  If 
too  few  angles  are  constrained,  the  results  are  generally  poor  since 
the  response  is  controlled  in  too  few  places,  (see  Appendix  C).  Simi¬ 
larly,  if  the  constrained  angles  are  concentrated  in  specific  areas, 
the  response  at  other  points  will  tend  to  be  quite  poor. 


It  has  been  found  that  a  loose  set  of  constraints,  such  as  incre¬ 
ment  angles  of  A0  =  5°  and  A«J»  =  20°  give  results  nearly  as  good  as 
those  for  a  tighter  set  of  constraints,  such  as  A0  =  2°  and  A4>  =5°, 
but  with  a  40%  to  75%  reduction  in  run  time.  The  range  of  reduction  in 
run  time  is  influenced  mainly  by  the  number  of  new  starts.  These  con¬ 
straints  are  valid  for  element  spacing,  d,  equal  to  V 2. 

It  is  obvious  that  with  tighter  constraints,  a  better  solution 
will  most  likely  be  found,  but  the  trade  off  is  an  unavoidable  increase 
in  run  time. 

4.4.3  MAXITE 

MAXITE  controls  the  maximum  number  of  iterations  of  the  Permuta¬ 
tion  Search  algorithm.  It  influences  both  the  amount  of  computer  time 
and  the  quality  of  the  final  results.  The  larger  MAXITE,  the  longer 
the  run  time.  For  smaller  arrays,  3  by  3  and  6  by  6,  20  iterations  was 
the  maximum  used  in  this  study.  Ten  iterations  was  the  maximum  used 
for  larger  arrays,  7  by  7  up  to  9  by  9.  Fewer  iterations  are  used  for 
larger  arrays  because  of  the  longer  run  time  associated  with  them. 
MAXITE  is  the  other  parameter  which  can  stop  program  execution.  The 
quality  of  result  is  affected  by  MAXITE  through  its  association  with 
the  effects  of  NCOUNT. 

4.4.4  NCOUNT 

NCOUNT  controls  the  maximum  number  of  new  starts  that  "PERMSER" 
will  use  during  its  search.  It  may  seem  that  setting  NCOUNT  very  high 
ensures  that  many  neighborhoods  will  be  searched,  thereby  improving  the 
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results.  In  practice,  this  is  not  the  case.  For  NCOUNT  values  greater 
than  5,  the  amount  of  improvement  for  each  new  start  decreases  with 
each  successive  new  start.  Thus  enters  the  problem  of  diminishing 
returns ! 

This  may  seem  to  imply  that  the  current  solution  is  near  the  glob¬ 
al  optimum.  Unfortunately,  this  is  not  usually  so.  If  two  arrays  are 
designed  from  the  same  selection  pool  of  simulated  elements,  the  re¬ 
sults  can  vary  greatly.  The  elements  selected  and  their  array  loca¬ 
tions  will  be  entirely  different.  If  both  attempts  are  random  starts, 
and  the  search  covers  the  same  number  of  neighborhoods,  the  resulting 
arrays  can  have  good  response  patterns,  or  bad  response  patterns. 
Therefore,  diminishing  returns  do  not  necessarily  herald  proximity  to  a 
global  optimum.  An  accurate  interpretation  is  that  the  quality  of  the 
locally  optimum  solution  depends  on  the  amount  of  area  in  the  solution 
space  that  has  been  searched. 

The  wide  variation  in  the  quality  of  responses  obtained  does  not 
discourage  the  use  of  Permutation  Search.  The  ratio  of  good  responses 
to  bad  responses  for  a  given  selection  pool  is  so  great  that  with  only 
a  few  runs  of  "PERMSER”,  an  acceptable  result  will  be  found. 


4.4.5  OTHER  INFLUENCES 


As  mentioned  earlier,  there  are  parameters  which  result  from  the 
problem  type  which  can  influence  "PERMSER” •  The  number  of  elements  in 
the  selection  pool  directly  influences  the  run  time,  and  to  a  lesser 
extent,  the  quality  of  results.  An  increase  in  the  number  of  elements 
will  increase  the  run  time  and  perhaps  improve  the  results  due  to  a 
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greater  number  of  permutations  available.  The  size  of  the  array  influ¬ 


ences  the  run  time  in  a  manner  similar  to  the  size  of  the  selection 


pool. 


4.4.6  THE  STARTING  POINT 


Several  tests  were  done  to  determine  the  effects  of  different 


starting  points.  The  different  types  of  starting  points  used  included 


random  arrangements,  sequential  arrangements,  and  group  arrangements  of 


elements  in  the  permutation  array.  Sequential  arrangements  are  those 


in  which  elements  adjacent  to  each  other  in  the  transducer  array  are 


adjacent  to  each  other  in  the  permutation  array.  Group  arrangements 


are  those  in  which  elements  that  occupy  symmetric  array  positions  are 


adjacent  in  the  permutation  array.  Many  variations  on  these  different 


arrangements  were  also  tried.  Results  indicate  that  the  best  array  re¬ 


sponse  patterns  are  obtained  from  randomly  arranged  initial  permuta¬ 


tions.  Those  arrays  which  result  from  group  arrangements  are  not  much 


poorer  than  those  from  the  random  arrangements.  They  differ  most  sig¬ 


nificantly  in  run  time  and  efficiency.  Efficiency  is  defined  as  the 


average  amount  of  improvement  per  base  permutation.  This  does  not  mean 


that  the  quality  of  the  initial  base  permutation  has  no  effect  on  the 


outcome.  It  means  the  arrangement  of  elements  in  the  permutation  array 


does  not  affect  the  objective  function  of  the  initial  base  permutation. 


The  most  significant  and  obvious  result  is  that  the  best  results 


are  obtained  from  those  initial  permutations  with  the  best  initial  re¬ 


sponse  patterns.  Therefore,  if  permutations  with  good  response  pat¬ 


terns  can  be  easily  found,  better  results  can  be  obtained. 


4.4.7  FLCW  CHART  OF  "PERMSER 
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The  following  is  a  discussion  of  the  operation  of  "PERMSER".  A 
flow  chart  for  "PERMSER"  is  shown  in  Figure  4.1.  The  element  para¬ 
meters  and  Permutation  Search  control  parameters  are  read  in  subroutine 
DATAINN.  Control  passes  from  DATAINN  to  NEARIN,  where  the  element 
parameters  are  written  into  an  output  file.  This  file  will  be  used  by 
"NEARPRES",  (see  Appendix  F),  to  compute  and  print  the  response  pattern 
for  the  initial  permutation.  The  value  of  the  objective  function  is 
determined  in  ZVALUE  and  is  used  as  the  initial  reference  point.  The 
interchanging  of  elements  occurs  in  PERMUTE.  Several  evaluations  are 
also  made  in  PERMUTE,  but  for  simplicity,  are  depicted  as  occuring  out¬ 
side  of  PERMUTE.  EVAL  determines  whether  the  trial  permutation  im¬ 
proves  the  objective  function.  To  do  this,  it  calls  ZVALUE.  ZVALUE 
calls  YVALUE  and  INITIALIZE  to  determine  the  response  of  each  trial 
permutation  at  the  constrained  angles.  INITIALIZE  determines  the  re¬ 
sponse  at  e  =  0°,  and  <}>  =  0°,  i.e.  ,  the  main  beam  response.  YVALUE  de¬ 
termines  the  response  at  every  constrained  angle  and  subracts  it  from 
the  main  beam  response  to  determine  how  many  dB  down  the  response  is  at 
each  constraint. 
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At  this  point,  there  are  two  possibilities.  Either,  the  trial 
permutation  becomes  a  new  base  permutation,  or  the  trial  permutation  is 
unsuccessful.  PERMUTE  will  continue  exchanges  of  elements  either  on 
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the  new  base  permutation,  or  on  the  previous  base  permutation.  When  a 
complete  cycle  of  exchanges  does  not  present  an  improvement,  NEWSTART 
is  called.  If  the  maximum  number  of  new  starts  has  not  been  exceeded, 
a  new  start  is  initiated.  If  the  maximum  number  has  been  exceeded,  the 
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Flowchart  for  "PERMSER 


search  process  stops.  Additionally,  after  a  complete  cycle  of  unsuc¬ 
cessful  exchanges  has  been  made,  PERMUTE  determines  if  MAXITE  has  been 
exceeded.  If  it  has,  the  search  process  is  terminated. 

Once  the  search  process  is  terminated,  the  results  are  printed 
out.  The  parameters  of  the  completed  array  are  written  into  a  file  to 
be  used  by  "NEARPRES”  to  print  a  response  pattern. 


CHAPTER  FIVE 


THE  NEO-KENDIG  SELECTION  SCHEME 

5.1  INTRODUCTION 

As  mentioned  in  Chapter  Four,  if  an  initial  base  permutation 
yields  a  good  response  pattern,  Permutation  Search  will  find  a  permu¬ 
tation  with  a  better  response  pattern.  The  problem,  then,  lies  in 
finding  a  good  starting  point. 

Since  the  Kendig  Scatter  Diagram  Method  has  been  used  success¬ 
fully  to  design  arrays  with  reduced  response  pattern  errors,  it  is  a 
logical  good  starting  point.  This  chapter  presents  a  new  selection 
scheme  based  on  the  Kendig  method.  Since  the  derivation  in  this  chap¬ 
ter  follow'-  Kendig's  derivation,  and  presents  the  same  idea  of  grouping 
elements  to  cancel  their  errors,  this  new  method  is  called  the  Neo- 
Kendig  Selection  Scheme.  However,  this  new  scheme  avoids  the  errors 
inherent  in  the  Kendig  method  by  using  a  three-dimensional  model  of  ar¬ 
ray  performance  instead  of  the  two-dimensional  model  used  previously. 

In  his  report  [1],  Kendig  presents  a  technique  to  minimize  the  im¬ 
pact  of  the  individual  transducer  element  errors  on  the  entire  trans¬ 
ducer  array  response.  In  essence,  the  Kendig  Scatter  Diagram  Method 
involves  the  paling  of  elements  which  will  occupy  symmetric  positions 
along  a  line  through  the  center  of  the  array  in  such  a  way  as  to  par¬ 
tially  cancel  their  respective  amplitude  and  phase  errors.  This  ch- 
nique  was  implemented  in  a  BASIC  routine  by  Schafer  [2].  It  represent¬ 
ed  an  improvement  over  the  more  tedious  hand  selection  version  previ- 


The  Neo-Kendig  method  represents  a  substantial  improvement  over 


the  existing  versions  of  the  Kendig  Scatter  Diagram  Method.  It  in¬ 
volves  the  simultaneous  selection  of  four  elements  to  occupy  symmetric 
array  positions,  instead  cf  pairs  of  elements  as  previously  done. 

5.2  DERIVATION 

In  this  derivation,  the  transducers  are  treated  as  point  elements 
located  at  the  center  of  the  transducer  faces.  The  contribution  of  an 
element,  m,  located  at  a  position  (Xm,Ym)  in  a  planar  array,  see  Figure 
5.1,  to  the  array  response  is 

Qm  =  V  exp(i(k*Um  +  am)}  ,  (5-1) 

where  Am  is  the  product  of  the  shading  for  position  m  and  the  amplitude 
response  of  the  transducer  element  in  that  position,  and  am  is  the  sum 
of  the  phase  shading  for  position  m  and  the  element's  phase  response. 
Um  is  a  function  of  the  x  and  y  coordinates  for  position  m  and  the  an¬ 
gles  of  incidence  6  and  4> >  viz., 

Um  =  Xm*cos<t>  sin9  +  Ym*sinb  sin9 

or  (5-2) 

Um  =  Vucos  +  Ym’usin  where  ucos  =  cos<f>  sin0 

usin  =  sin«j)  sin0  . 

Expanding  the  exponential  into  sin  and  cos  form  yields 

exp{i(k*Um  +  am)}  =  cos(k’Um)cos(am)  -  sin(k*Um)sin(am)  (5-3) 

+  i[sin(k*Um)cos(am)  +  cos(k*Um)sin(am) ] . 


Using  Equation  5-2,  the  sin  and  cos  terras  of  Equation  5-3  can  be 


further  expanded  into 


cos(k*Um)  =  cos(k*Xm*ucos  +  k*Ym*usin) 

=  cos (k'X^'ucos )cos (k* Ym*usin) 

-  sin(k*Xm*ucos )sin(k* Ym*usin) 


( 5-4) 


sin(k*Um)  =  sin(k*Xm*ucos  +  k*Ym‘usin) 

=  sin(k*Xm*ucos )cos (k* Ym*usin) 

+  cos(k,Xm*ucos)sin(k*Ym*usin). 


(5-5) 


Consider  four  positions  located  symmetrically  about  the  center  of 
an  even  array,  (Figure  5.2).  The  contribution  of  these  four  elements 


A| *exp(i(k*Uj  +  a])}  +  A2'exp{i(k*U2  +  a2)}  (5-6) 

+  A3*expti(k*U3  +  83))  +  A4’exp(i(k'U4  +  84)}, 


where 


Uj  =  X^’ucos  +  Yj^usin  U3  =  -Xj’ucos  -  Yjusin 
u2  =  “Xj'ucos  +  Yjusin  U4  =  X^*ucos  -  Yjusin 


(5-7) 


Note  that  Uj  =  -U3  and  U2  =  -U4  .  It  is  assumed  there  are  no  position¬ 
al  errors  when  the  elements  are  placed  in  the  array. 

Breaking  the  response  up  into  real  and  imaginary  components,  the 
real  part  can  now  be  expressed  as 


Re{Q4}  = 


[cos(k*Ui )cosa[  -  sin(k*Uj )sinaj ]  + 
A2(cos(k*U2)cosa2  ~  sin(k’U2)sina2]  + 
A3(cos(k*U3)cosa3  -  sin(k*U3)sina3]  + 
A4[cos(k*U4)cosa4  -  sin(k*U4)sina4) . 


(5-8) 
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For  small  phase,  l.e.,  a  «  I,  (in  radians),  the  following  approxima¬ 
tions  can  be  made, 

cos  a-L  =  1  and  sin  a^  =  a^  .  (5-9) 

Substituting  Equations  5-7  and  5-9  into  Equation  5-8,  the  real  and 
imaginary  parts  of  the  contribution  can  be  written  as  shown  in  Equa¬ 
tions  5-10  and  5-11,  (Figure  5.3).  The  terms  in  square  brackets  will 
be  referred  to  as  the  Kendig  Terms  and  distinguished  by  the  number  ap¬ 
pearing  above  each  term.  The  sin  and  cos  products  will  be  referred  to 
as  the  Kendig  Coefficients. 

This  preceding  derivation  is  valid  for  all  four-fold  symmetric 
even  arrays.  The  derivation  for  odd  arrays  does  not  differ  substan¬ 
tially  and  produces  an  equivalent  result.  It  is  presented  in  Appendix 

B. 

5.3  THE  KENDIG  TERMS 

In  Equations  5-10  and  5-11,  terras  3,  4,  7,  and  8  include  both  the 
phase  and  amplitude  responses  of  the  elements,  while  terms  1,  2,  5,  and 
6  represent  contributions  from  amplitude  only.  If  all  the  elements 
have  ideal  response  characteristics,  e.g.,  =  1.0  and  =  0.0, 

only  one  term,  term  1,  would  be  non-zero.  Thus,  terms  2  through  8  rep¬ 
resent  the  error  Introduced  into  the  array  response.  Minimizing  the 
total  error  entails  minimizing  all  seven  of  these  terms. 


A  look  at  terms  3  and  4  reveals  a  dilemma.  Minimizing  terra  3 


means 


or 


A3a3  +  A4a4  =  A^a^  +  A4a4 
A3a3  ~  Alal  =  A2a2  "  A4a4  » 


while  minimizing  term  4  means 


or 


A3a3  +  A2a2  =  A^a^  + 

A3a3  -  Alal  =  a4^4  -  A2a2  • 


(5-12) 


(5-13) 


Excluding  the  uncommon  case  of  A2a2  =  k^a^,  minimizing  term  3  will  nec¬ 
essarily  prevent  an  equal  reduction  of  the  value  of  term  4.  This  rela¬ 
tionship  between  terms  3  and  4  exists  between  all  combinations  of  terms 
3,  4  and  8,  as  well  as  combinations  of  terms  2,  5,  and  6,  (see  Table 
5.1). 

Consider  the  Kendig  Coefficients.  For  any  array,  x  and  y  are 
known  for  every  position.  Since  ucos  and  usin  are  functions  of  the  in¬ 
cident  angles,  the  Kendig  Coefficients  can  be  determined  for  specific 
incident  angles.  If  elements  are  then  chosen  in  such  a  way  as  to  mini¬ 
mize  the  sum  at  each  symmetric  position  of  the  product  of  the  terms  and 
coefficients,  the  error  at  this  angle  is  reduced.  However,  at  a  dif¬ 
ferent  incident  angle  which  is  not  nearby  the  first,  this  minimization 
no  longer  necessarily  holds  because  the  incident  angles  result  in  dif¬ 
ferent  Kendig  Coefficient  values.  For  a  more  detailed  discussion  of 
the  implications  of  this,  see  Appendix  C. 


The  solution  to  reducing  the  error  impact  lies  in  reducing  each 
terra  independently  of  the  others,  (in  this  case,  independently  means 
without  regard  to  the  interaction  of  the  Kendig  Coefficients).  If  each 
terra  is  as  small  as  possible,  the  product  with  the  coefficients  will 
also  be  small.  Thus,  the  error  impact  will  be  reduced  over  all  angles. 

5.4  THE  SELECTION  ROUTINE 

As  just  discussed  in  the  previous  section,  reducing  the  total  er¬ 
ror  introduced  into  the  array  response  involves  minimizing  the  Kendig 
Terms  individually.  In  order  to  determine  which  terms  to  minimize  and 
in  what  order,  the  Kendig  terras  of  Equations  5-10  and  5-11  were  pro¬ 
grammed  into  a  BASIC  routine  on  an  HP9825B  calculator.  This  routine 
requires,  as  an  input,  the  radiation  loading  for  each  symmetric  array 
position.  Using  these  and  the  equivalent  circuit  parameters  of  each 
element,  it  computes  the  amplitude  and  phase  response  for  every  element 
in  each  symmetric  position.  Using  one  of  the  Kendig  terms  as  a  criter¬ 
ion,  the  routine  searches  through  the  set  of  responses  for  a  chosen  po¬ 
sition,  and  finds  the  group  of  four  that  gives  the  lowest  value  for  the 
specified  term.  Since  terras  3,  4,  7,  and  8  represent  the  error  due  to 
both  the  phase  and  amplitude  responses  of  an  element,  combinations  of 
these  terras  were  the  Kendig  terras  used.  Although  terms  2,  5,  and  6  ap¬ 
pear  to  depend  on  amplitude  only,  the  phase  responses  are  still  repre¬ 
sented  in  these  terms  through  the  approximation  cos  =  1.  These 
terras  are  not  used  because  the  phase  influence  is  much  less. 
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5.4.1  TEST  RESULTS 


Many  different  tests  of  the  selection  process  were  done.  Using 
the  results  of  section  2.1.2,  several  selection  pools  were  generated  on 
the  computer.  Different  combinations  of  terms  3,  4,  7,  and  8  were  used 
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to  test  the  selection  process.  These  selection  pools  were  used  in  var¬ 
ious  ways  to  design  many  arrays.  Combinations  of  one,  two  ,  three  or 
four  term(s)  were  used  for  the  selection  criteria. 
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It  has  been  found  that  use  of  any  two  terms  for  the  selection  cri¬ 
teria  gives  the  best  results.  It  is  not  surprising  that  minimizing  two 
terms  would  give  better  results  than  minimizing  only  one.  However,  it 
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is  surprising  that  minimizing  three  or  four  terms  does  not  present  add¬ 
itional  improvement.  There  are  no  apparent  causes  for  this  phenomenon. 
For  a  collection  of  selection  pools,  there  is  no  particular  com¬ 


bination  of  terms  3,  4,  and  8  that  will  give  the  best  response  for  ev¬ 
ery  selection  pool.  Any  combination  of  terms  will  choose  the  same  four 


elements  for  each  symmetric  position,  but  will  arrange  them  differently 


within  the  array.  These  different  arrangements,  or  permutations,  can 
sometimes  result  in  vastly  different  responses,  and  other  times  give 
comparatively  equivalent  responses.  In  addition,  if  two  arrays  are 
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created  from  a  large  selection  pool,  the  combination  of  terms  which  re¬ 
sults  in  a  good  first  array  does  not  necessarily  produce  a  good  second 
array. 


Term  7  has  been  excluded  from  consideration  because  it  does  not 
have  the  same  form  as  terras  3,  4,  and  8;  it  offers  no  information  on 
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the  arrangement  of  elements  in  the  groups  selected.  Additionally,  for 
small  amplitude  and  phase  variances,  the  range  of  values  for  term  7  is 
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not  as  great  as  for  the  other  terms.  Greater  reduction  of  error  impact 
can  be  achieved  by  striving  for  the  minimum  value  of  the  other  terms 
than  for  the  minimum  value  of  term  7. 

5.5  THE  NEO-KENDIG  COMPUTER  ROUTINE,  "HKV" 

The  following  is  a  brief  description  of  the  actual  selection  pro¬ 
cess.  For  an  even  array,  the  radiation  loading  for  the  symmetric  posi¬ 
tion  closest  to  the  array  center  is  used  to  determine  the  amplitude  and 
phase  responses  of  all  the  elements  in  the  selection  pool.  For  these, 
the  means  and  variances  are  found.  Elements  that  are  farthest  from  a 
reference  point  are  temporarily  removed  and  new  means  and  variances  are 
found.  Again,  the  elements  farthest  from  the  means  are  removed.  This 
process  is  repeated  until  there  are  between  10  and  14  elements  remain¬ 
ing  in  the  pool. 

Several  different  methods  of  searching  the  selection  pool  for  the 
10  to  14  initial  elements  were  also  tried.  These  included  using  the 
initial  means  of  the  entire  selection  pool  as  the  reference  point  for 
the  entire  array,  or  using  no  reference  point.  It  was  determined  that 
use  of  the  first  four  chosen  elements  as  a  reference  point  for  all  sub¬ 
sequent  positions  produces  the  best  results.  The  reference  point  to 
choose  the  first  four  is  the  moving  mean  just  described. 

The  value  of  the  first  Kendig  terra  is  then  determined  for  every 
possible  permutation  of  four  elements,  a  total  of  24024  trials  for  14 
elements.  The  lowest  40  of  these  are  kept.  The  value  of  the  second 
term  is  then  ascertained  for  these  40  best  permutations.  The  permuta¬ 
tion  that  produces  the  minimum  response  for  the  second  term  is  the 


permutation  for  the  first  position.  Next,  the  chosen  four  elements  are 
permanently  removed  from  the  selection  pool  and  those  temporarily  re¬ 
moved  are  restored. 

The  means  and  variances  of  these  four  elements  are  determined  and 
used  as  a  reference  point  for  all  subsequent  selections.  The  selection 
process  for  every  successive  position  follows  the  same  procedure  except 
that  the  elements  removed  are  those  farthest  from  the  means  of  four 
elements  first  chosen.  Consequently,  it  is  not  necessary  to  recalcu¬ 
late  the  means  and  variances  of  the  elements  in  the  selection  pool, 
once  the  first  symmetric  position  is  filled. 

The  routine  for  an  odd  array  is  different  in  one  respect.  Since 
the  array  has  one  unsymmetric  position  at  the  center  of  the  array,  the 
element  nearest  to  means  of  the  group  of  10  to  14  elements  is  chosen 
and  placed  in  the  center  position.  The  amplitude  and  phase  response  of 
this  element  is  used  as  the  reference  point  for  all  subsequent 
selections . 

5.5.1  FLOW  DIAGRAM  FOR  "HKV" 

The  Neo-Kendig  method  has  been  implemented  into  a  FORTRAN  computer 
routine  called  "HKV”,  (a  listing  is  in  Appendix  E).  This  routine 
allows  any  combination  of  Kendig  Terras  3,  4,  and  8  to  be  used  as  the 
selection  criteria.  A  flow  diagram  is  shown  in  Figure  5.4. 

The  element  circuit  parameters,  and  selection  control  parameters 
are  read  in  by  the  routine  DATAININ.  The  element  responses  for  the 
first  symmetric  array  position,  (always  closest  to  the  center),  are 
computed  in  RESPONSE.  The  mean  amplitude  and  phase  responses  for  the 
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selection  pool  are  computed  and  stored  for  reference  in  STAT.  Those 
elements  which  are  farther  than  a  specified  tolerance  from  the  selec¬ 
tion  pool  means  are  permanently  removed  from  the  selection  pool  in 
DELETE. 

Routines  REMOVE  and  OPTIMO  temporarily  remove  unsuitable  elements 
from  the  selection  pool  and  search  for  those  permutations  which  give 
minimal  values  for  the  Kendig  Terms  used  as  the  search  criteria. 
OPTIMO  calls  USET  to  determine  the  values  of  the  Kendig  Terms  being  an¬ 
alyzed.  Once  such  a  permutation  is  found,  the  elements  are  assigned  to 
their  respective  array  positions,  and  permanently  removed  from  the  se¬ 
lection  pool.  If  the  array  position  just  filled  is  the  centermost  po¬ 
sition,  the  chosen  elements'  means  are  computed  and  stored  for  refer¬ 
ence  in  SETK.  In  PREPNSET,  those  elements  temporarily  removed  are  now 
returned  to  the  selection  pool  to  be  considered  for  the  next  array  po¬ 
sition.  The  process  loops  back  to  RESPONSE  to  calculate  the  element 
responses  for  the  next  array  position,  and  the  selection  process  re¬ 
peats.  Upon  completion  of  the  selection  process,  WRITEIT  prints  a  map 
of  the  completed  array,  showing  where  each  element  is  situated.  Rou¬ 
tines  NEARIN1  and  PERMSERIN  are  used  to  write  the  parameters  of  the 
completed  array  into  files  for  use  by  "NEARPRES”  and  "PERMSER”. 


CHAPTER  SIX 


THE  SELEG..ON  PROCESS 


6.1  THE  COMBINATION 

As  mentioned  in  Chapter  Four,  Permutation  Search  will  look  for 
better  permutations  in  the  neighborhood  of  an  initial  permutation. 
However,  without  a  good  starting  point,  the  results  are  not  always  ac¬ 
ceptable.  With  the  Neo-Kendig  method,  good  results  can  be  found,  but 
not  guaranteed.  Several  different  combinations  of  two  terms  must  be 
tried  in  order  to  find  an  acceptable  result. 

Combining  the  Permutation  Search  and  Neo-Kendig  methods  results  in 
a  cancellation  of  their  respective  disadvantages.  If  the  solutions  ob¬ 
tained  from  the  Neo-Kendig  method  are  used  as  initial  permutations  for 
"PERMSER",  the  results  are  superior  to  those  obtained  by  either  method 
alone.  Any  pairing  of  terms  3,  4,  and  8  produces  results  which  are  ac¬ 
ceptable  initial  permutations  for  "PERMSER”.  As  discussed  in  Chapter 
Four,  a  random  arrangement  of  elements  in  the  permutation  array  gives 
better  results  than  placing  elements  occupying  symmetric  positions  ad¬ 
jacent  to  each  other.  The  results  of  the  methods  described  in  Chapters 
Four  and  Five  are  presented  in  this  chapter  to  facilitate  comparison. 
Figure  6.1  shows  one  half  of  the  response  pattern  for  a  6  by  6  array 
of  ideal  elements  and  a  desired  side  lobe  level  of  -40  dB.  The  shading 
values  are  from  the  Do lph-Chebyshev  method  for  line  arrays  [12]  and  the 
second  product  theorem  [13].  This  map  represents  the  response  of  the 


array  to  a  signal  at  any  incident  angle,  9,  between  0°  and  -90°  and 
roll  plane,  4>  »  between  0°  and  180°.  The  reference  level  is  the  re¬ 
sponse  of  the  array  to  the  same  signal  at  9  =  0°  and  <{>  =  0°  incidence. 
The  dB  levels  represesent  negative  decibels.  A  colon,  (:),  following 
a  number  indicates  a  phase  greater  than  90°  different  from  the  response 
at  the  peak  of  the  main  beam,  referred  to  as  a  phase  change.  Figure 
6.2  shows  the  same  portion  of  the  response  pattern  for  a  random  ar¬ 
rangement  of  non-ideal  elements.  The  selection  pool  was  the  set  of  66 
transducers  mentioned  in  Chapter  Two.  The  regions  of  the  response  pat¬ 
tern  that  meet  the  40  dB  level  requirement  are  outlined.  Note  how  the 
locations  of  the  phase  changes  have  been  altered  significantly.  Pre¬ 
sented  in  Figure  6.3  is  the  result  obtained  after  17  iterations  of  the 
Permutation  Search  algorithm  for  9  constrained  every  5°  and  <$>  con¬ 
strained  every  20°  in  the  side-lobe  region.  This  result  was  obtained 
from  a  random  initial  permutation.  There  are  only  a  few  locations 
which  do  not  meet  the  desired  40  dB  level.  In  addition  to  the  amount 
of  area  that  meets  the  40  dB  requirement,  a  good  indication  of  the 
quality  of  the  result  is  that  the  locations  of  the  phase  changes  ap¬ 
proximate  those  of  the  ideal  case. 

The  same  selection  pool  that  produced  the  response  just  mentioned 
was  also  used  by  the  Neo-Kendig  Selection  Scheme.  Terms  3  and  4  were 
the  selection  criteria.  The  results  are  shown  in  Figure  6.4.  Again, 
the  locations  of  the  phase  changes  approximate  those  of  the  ideal  ar¬ 
ray.  In  this  case,  the  worst  response  Is  38  dB,  only  a  2  dB  error. 


m 
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Finally,  Figure  6.5  shows  the  response  pattern  obtained  by  using 
the  results  of  the  Neo-Kendig  method  as  the  initial  starting  point  for 
"PERMSER".  The  only  deviations  from  the  desired  40  dB  level  are  1  dB 
in  magnitude  and  not  concentrated  in  one  area.  This  result  was  ob¬ 
tained  after  12  iterations,  during  which  two  new  starts  were  dene  with 
a  run  time  of  31  epu  minutes  on  a  VAX  11/782. 

Similar  results  for  an  8  by  8  array  with  three  elements  removed 
from  each  corner  are  shown  in  Figures  6.6  -  6.9.  The  selection  pool 
for  this  array  contained  140  elements  and  was  generated  as  described  in 
Chapter  Two.  The  shading  values  for  the  desired  side  lobe  levels  of 
-40  dB  were  obtained  by  Wilson  [10].  Figure  6.6  is  the  response  for 
ideal  elements.  Figure  6.7  is  for  a  random  arrangement  of  elements. 
Figure  6.8  is  the  response  pattern  obtained  by  the  Neo-Kendig  method. 
Figure  6.9  is  the  final  result,  obtained  after  11  iterations  of  the 
Permutation  Search  algorithm,  three  new  starts,  and  80  epu  minutes. 

The  maximum  tolerance  on  the  elements  was  9%  in  amplitude,  and  11° 
in  phase.  However,  the  tolerances  on  elements  in  the  selection  pool  do 
not  indicate  which  elements  will  be  selected.  The  Neo-Kendig  method 
will  use  only  those  elements  which  are  close  to  the  reference  means. 
Therefore,  those  elements  farthest  from  the  means  will  not  be  selected. 
But,  "PERMSER"  disregards  an  element's  distance  from  a  reference  mean. 
Any  element  may  enter  the  array  if  its  presence  produces  better  re¬ 
sults.  It  is  possible,  and  has  happened,  that  elements  with  large  de¬ 
viations  will  be  placed  in  an  array  alongside  elements  with  small 
deviations . 
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CHAPTER  SEVEN 


SUMMARY  AND  CONCLUSIONS 


7 . 1  SUMMARY 

The  goal  of  this  thesis  was  to  develop  an  element  selection  tech¬ 
nique  which  would  place  elements  in  an  array  in  a  manner  that  minimizes 
the  impact  of  individual  transducer  element  errors  on  the  array's  di¬ 
rectional  beam  pattern.  The  combination  of  this  technique  and  the  ac¬ 
curate  measurement  system  developed  by  Schafer  [2]  is  a  powerful  tool 
which  can  reduce  costs  and  production  time  during  the  development  of  a 
transducer  array. 

Chapter  Two  discusses  various  aspects  of  the  lumped-parameter 
equivalent  circuit.  The  correlations  between  the  measured  parameters 
have  been  established.  A  preliminary  set  of  tolerances  on  the  measured 
circuit  parameters  has  been  established,  based  on  a  measured  set  of 
transducers.  These  tolerances  indicate  which  circuit  parameters  have 
the  most  influence  on  an  element's  response. 

Chapter  Three  discusses  the  mutual  interaction  of  transducer  ele¬ 
ments  in  an  array.  The  effect  of  incidence  angle  is  discussed,  as  well 
as  the  effects  of  element  tolerances  on  the  mutual  interaction. 

The  application  of  Permutation  Search  is  described  in  Chapter 
Four.  The  control  parameters  and  the  type  of  results  that  can  be  ex¬ 
pected  are  discussed.  The  Neo-Kendig  selection  technique  is  derived  in 
Chapter  Five.  It  is  an  effective  element  placement  technique  by 


itself.  When  it  is  used  as  a  starting  point  for  the  Permutation  Search 
algorithm,  the  results  are  spectacular.  This  combination  has  designed 
arrays,  6  by  6  and  8  by  8,  using  elements  with  maximum  tolerances  of 
*9%  in  amplitude  and  ±11°  in  phase,  that  meet  a  required  40  dB  side 
lobe  level  in  all  but  a  very  few  locations,  and  fail  at  those  points  by 
less  than  2  dB. 

7.2  FUTURE  CONSIDERATIONS 

Permutation  Search  is  the  most  time  consuming  part  of  the  entire 
placement  process.  A  reduction  of  the  run  time,  without  a  concurrent 
loss  of  effectiveness  or  a  reduction  of  result  quality,  would  represent 
an  additional  improvement  of  the  process. 

The  derivation  of  the  Neo-Kendig  technique  has  shed  light  on  a 
mechanism  that  plagues  the  design  of  shading  values.  The  use  of  the 
Kendig  Coefficients  for  a  given  array  may  reduce  the  time  necessary  for 
developing  shading  values  by  eliminating  some  computations.  Addition¬ 
ally,  if  the  shading  values  are  considered  variables,  and  the  Kendig 
Coefficients  constants,  the  problem  of  achieving  a  desired  beam  pattern 
becomes  a  set  of  linear  equations  which  may  be  solvable  by  linear 
programming. 

The  selective  use  of  weighting  may  improve  the  effectiveness  of 
Permutation  Search.  By  assigning  greater  weighting  values  at  regions 
in  the  response  plane  where  it  is  more  crucial  that  the  specified  side 
lobe  level  be  attained,  greater  control  over  these  areas  may  be 


achieved. 
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A  similar  result  might  be  achieved  through  the  use  of  Fuzzy  Math, 
[9].  Assigning  membership  functions  to  different  regions  of  the  re¬ 
sponse  pattern  would  provide  more  information  on  the  acceptabilty  of 
the  current  base  permutation  than  is  provided  by  the  value  of  the  ob¬ 
jective  function.  For  example,  Fuzzy  Math  would  allow  the  case  where 
achieving  a  response  better  than  40  dB  in  a  certain  region  would  out¬ 
weigh  the  negative  effects  of  not  achieving  40  dB  in  another  region. 

Different  methods  of  error  analysis  may  improve  both  the  quality 
of  results  and  run  time.  One  example  would  be  to  minimize  the  maximum 
error  rather  than  the  total  error,  as  is  the  case  in  this  work. 

The  effects  of  incident  angle  on  the  interactions  between  elements 
was  mentioned  in  Chapter  Two.  If  the  assumption  of  nominal  phase  ef¬ 
fects  could  be  validated,  this  selection  process  would  stand  on  firmer 
ground.  For  every  incident  angle,  the  mutual  interactions  present  N 
equations  and  N  unknowns,  where  N  is  the  array  size.  Through  matrix 
manipulation,  it  is  possible  to  solve  for  all  the  interaction  coeffi¬ 
cients,  and  confirm  or  refute  that  assumption. 

The  testing  of  more  sets  of  transducer  elements  would  enable  the 
establishment  of  more  reliable  equivalent  circuit  parameter  tolerances. 
This,  in  turn,  would  allow  more  uniformity  in  element  performance 
characteristics  to  be  achieved.  With  increased  uniformity  comes  in¬ 
creased  array  performance. 


APPENDIX  A 


ARRAY  PARAMETERS  AND  TERMINOLOGY 
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The  arrays  considered  in  this  study  are  four-fold  symmetric,  (or 
quadrature  symmetry),  even  and  odd  arrays.  An  odd  array  is  one  which 
has  an  element  at  the  center  of  the  array,  at  X^  =  0.0,  Yj  =0.0.  An 
even  array  has  no  center  element.  There  are  several  parameters  used  to 
describe  these  types  of  arrays.  Accompanying  this  description  are  ex¬ 
amples  of  an  even  and  an  odd  array.  In  order  to  implement  easily  the 
Neo-Kendig  Selection  Scheme  and  the  Permutation  Search  algorithm,  it 
was  necessary  to  create  an  array  numbering  system  which  reflects  the 
unique  symmetry  of  these  arrays. 

The  parameter  M  describes  the  number  of  positions  in  a  full  array, 
i.e.,  with  no  corner  elements  removed.  The  number  of  positions  in  one 
quadrant  is  described  by  IM.  For  even  arrays,  IM  =  M/4.  For  odd  ar¬ 
rays,  IM  =  (M  -  l)/4.  IM1  is  the  number  of  positions  in  one  quadrant 
that  will  be  occupied  by  a  transducer  element.  If  three  elements  are 
removed  from  each  corner,  then  IM1  =  IM  -  3. 

Each  position  in  the  array  has  a  unique  number  assigned  to  it. 
The  array  is  numbered  sequentially,  as  shown  in  Figures  A. 1  and  A. 2. 
For  the  purposes  of  assigning  shading  values,  each  position  in  the 
first  quadrant,  (and  the  center  element  in  an  odd  array),  is  assigned  a 
letter  representing  the  appropriate  shading  value.  Examples  of  shading 
values  are  given  for  desired  side  lobe  levels  of  -30  dB  and  -40  dB. 
These  are  obtained  from  the  Dolph-Chebyshev  [12]  shading  values  for  a 
line  array  and  the  second  product  theorem  [13]. 
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M  =  64 

10 

11 

12 

IM  =  16 

,1 

6 

7 

8 

IM1  =  13 

II 

elements  are 
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Side  lobe  level 


Position 

A 

B 

C 

D 

E 

F 

G 

H 

J 

K 


-30dB 

1.0000 

.8120 
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.0687 


-40dB 

1.0000 
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.0610 
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Figure  A.l 

Parameters  of  an  Even  Array,  8  by  8 
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Side  lobe  level 


Position 

A 

B 

C 

D 

E 

F 

G 

H 

J 

K 


-30dB 

1.0000 
.8738 
.5683 
.2642 
.  7635 
.4966 
.2309 
.  3230 
.1501 
.0698 


1.0000 
.8397 
.4793 
.1594 
.  7051 
.4025 
.1338 
.2297 
.0764 
.0254 


Figure  A. 2 

Parameters  of  an  Odd  Array,  7  by  7 
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The  term  "symmetric  array  position"  refers  to  those  elements,  one 


in  each  quadrant  which  have  the  relationship  described  in  equation  5-7, 


(even),  or  equation  B-l,  (odd).  In  the  even  array,  (this  holds  true 


for  the  odd  array  as  well),  the  symmetric  array  position  numbers  are 


arithmetically  separated  by  the  value  of  IM.  For  example,  for  the  8  by 


8  array  in  Figure  A. 1 ,  positions  2,  18,  34,  and  50  are  symmetric* 


Their  arithmetic  differences  are  16,  which  is  IM  for  an  8  by  8  array. 


All  symmetric  array  positions  have  the  same  shading  values. 
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APPENDIX  B 


THE  ADJUSTMENT  TO  NEO-KENDIG  FOR  ODD  ARRAYS 

There  are  two  types  of  symmetry  present  in  an  odd  array.  Consider 
the  first  two  quadrants  of  positions  in  an  odd  array. 


Figure  B.I 

First  Two  Quadrants  of  an  Odd  Array 

The  first  type  is  reflective  in  nature.  Position  6  is  symmetric  about 
the  y-axis  with  position  20. 

(x2o>y2o)  =  (-x6.V  (B_1> 

Those  positions  which  do  not  lie  along  the  axes  of  the  array  have  three 
other  symmetric  positions,  for  a  total  of  four.  However,  difficulties 
arise  with  this  type  of  symmetry  when  considering  position  16.  Since 
it  lies  along  one  of  the  axes,  (the  y-axis  in  this  case),  there  is  only 
one  other  position  symmetric  with  it. 
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The  other  type  of  symmetry  present  is  rotational.  Position  6  is 
symmetric  with  position  18  through  a  90°  rotation  about  the  origin. 
The  symmetry  is  such  that  the  x  and  y  coordinates  are  interchanged, 
viz.  , 


(Xi8»*18)  =  <>y6>x6)' 


(B-2) 


Note  that  the  arithmetic  relationship  between  these  two  positions  is 
the  same  as  between  two  symmetric  positions  in  an  even  array,  i.e., 
18  -  6  =  12,  which  is  IM  for  a  7  by  7  array.  Figure  B.2  is  the  odd 
array  equivalent  of  Figure  5.2. 

Since  the  Neo-Kendig  Selection  Scheme  takes  advantage  of  the 
unique  symmetry  of  an  even  array,  it  would  be  advantageous  to  extend 
the  selection  scheme  to  odd  arrays  with  as  little  change  as  possible. 
A  selection  scheme  for  odd  arrays  based  on  the  first  type  of  symmetry 
would  not  resemble  the  scheme  for  even  arrays  since  there  are  positions 
which  are  not  four-fold  symmetric.  The  second  type  of  odd  symmetry 
more  closely  resembles  the  even  symmetry  of  Chapter  Five.  It  is  this 
type  of  symmetry  which  is  used  for  odd  arrays. 

Paralleling  the  derivation  of  Chapter  Five,  equation  B-3  is  the 
odd  array  form  of  equation  5-7. 


=  Xj’ucos  +  Y^'usin  U3  =  -Xj'ucos  -  Y^’usin 

U2  =  -Y^'ucos  +  Xj*usin  U4  =  +Yj*ucos  -  Xj’usin 


(B-3) 


Therefore , 


UL  =  -U3 


u2  -  -U4 


(B-4) 


is  the  same  as  for  an  even  array. 
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The  difference  between  an  even  and  an  odd  array  shows  up  in  the 
odd  array  equivalent  of  equations  5-4  and  3-5,  which  represent  the  Ken- 
dig  Coefficients.  Since  the  coefficients  do  not  affect  the  selection 
process,  the  rest  of  the  derivation  is  not  presented. 

The  Neo-Kendig  Selection  Scheme  for  even  arrays,  outlined  in  Chap¬ 
ter  Five,  can  be  applied,  with  a  small  modification  to  deal  with  the 
center  element,  directly  to  odd  arrays. 

The  Neo-Kendig  technique  is  based  on  the  inherent  symmetry  of  the 
arrays  considered.  If  an  array  with  a  different  type  of  symmetry  is 
considered,  the  results  of  Chapter  Five  are  not  applicable.  However,  a 
similar  derivation,  reflecting  the  different  symmetry,  would  yield  a 
Neo-Kendig  technique  for  that  symmetry.  The  derivation  is  straight¬ 
forward  and  not  lengthy. 


APPENDIX  C 


KENDIG  COEFFICIENTS  AND  ARRAY  RESPONSE 

The  Dolph-Chebyshev  [12]  method  is  a  well  known  technique  used  to 
determine  the  shading  values  to  give  an  optimally  narrow  beam  width  for 
a  specified  uniform  side-lobe  level  for  line  arrays.  The  second  pro¬ 
duct  theorem  [13]  extends  this  result  to  planar  arrays.  However,  the 
shading  values  for  other  types  of  response  patterns  are  not  generated 
as  easily.  Computer  routines  which  use  linear  and  goal  programming 
techniques  can  be  used  to  generate  these  values  [10],  [11]. 

The  use  of  these  programs  has  revealed  a  curious  phenomenon.  Con¬ 
straining  a  certain  area  of  the  response  pattern  to  a  very  low  side- 
lobe  level  will  result  in  an  increase  in  side-lobe  level  at  another 
area  in  the  response  pattern,  resembling  that  which  happens  to  a  water 
balloon  when  pushed  in  at  one  point;  it  bulges  out  at  another  point. 
The  following  is  a  description  of  the  mechanism  for  this  phenomenon,  as 
enlightened  by  the  Neo-Kendig  coefficients. 

Without  loss  of  generality,  the  following  assumptions  can  be  made. 

1)  Only  ideal  elements  are  considered. 

2)  There  is  no  phase  shading. 

3)  k  =  1.0. 

Term  1  is  the  only  non-zero  Kendig  Term.  For  a  6  by  6  array  and  an  in¬ 
cident  angle  of  4>  =  75°  and  0  =  45°,  the  Kendig  Coefficients  for  each 


IUW  W*  IWmWMW  JJ'AWJV  W  WIW  g™  ^  ^ j- 


JT?’>’J<  V  V  V  V  V 


95 


symmetric  position  are  shown  in  Figure  C.l.  Each  of  the  coefficients 
for  the  symmetric  positions  is  multiplied  by  four  times  the  shading 
value  for  that  position,  (once  for  each  quadrant).  The  response  of  the 
array  is  twenty  times  the  logarithm,  base  10,  of  the  sum  of  these  pro¬ 
ducts  over  the  symmetric  positions  with  the  main  beam  response  as  the 
reference  level.  If  it  is  assumed,  for  example,  that  the  response  at 
this  incident  angle  is  very  low,  i.e.,  a  null,  the  sum  of  these  pro¬ 
ducts  over  all  the  symmetric  positions  is  small.  The  logarithm  of  a 
small  number,  <  1 ,  is  a  negative  number,  and  in  this  case,  represents 
the  dB  down  from  the  main  beam  response. 

At  a  different  incident  angle,  <*>  =  5°  and  0  =  60°,  which  is  not  in 
the  neighborhood  of  the  first,  the  Kendig  Coefficients  at  each  sym¬ 
metric  position  will  have  different  values,  (Figure  C.2).  Therefore, 
the  sum  of  the  products  at  the  second  angle  will  not  be  the  same  as  the 
sum  for  the  first  angle.  If  the  first  sura  was  very  small,  i.e.,  near 
zero,  this  second  sum  will  be,  in  most  cases  including  this  example. 


greater  than  the  first.  The  logarithm  of  this  larger  sum  will  result 
in  a  response  level  which  is  not  as  low  as  for  the  smaller  sum. 

This  explains  why  driving  the  response  to  a  very  low  level  in  one 
region  of  the  response  plane  may  degrade  the  response  in  another  re¬ 
gion. 
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THETA  =  60 
ustn  =  .07  55 


PHI  =5 
ucos=.8627 


Figure  C.2 


Kendig  Coefficients  for  0  =  60°, 
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APPENDIX  D 


LISTING  OF  PERMUTATION  SEARCH  COMPUTER  PROGRAM 


c**** ***************************************************** ************** 

PERMSER  * 

PERMUTATION  SEARCH  COMPUTER  ALGORITHM  * 

WRITTEN  BY  KING  W.  WIEMANN  * 

C*** ********************************** ********************************** 

COMMON/ COMMO 2/  Y 

COMMON/ COMMO 3/  RER, IMR,RERR 

COMMON/ COMM04/  K,M,AK 

COMMON/ COMMO 5/  IEV 

COMMON/ COMMI1/  NACH, NPRIOR, NOB J 

COMMON/ COMM 12/  ZBEST 

COMMON/COMM17/  MM,NPOOL,lJ 

DIMENSION  ZBEST( 10) , Y( 3000) , IMR(200) ,Z( 10) 

INTEGER  RER( 200) , RERR( 200) 

CALL  DATAINN 

C  DETERMINE  THE  INITIAL  ACHIEVEMENT  LEVEL. 

CALL  ZVALUE(Z) 

DO  10  KK=1 , NPRIOR 
ZBEST(KK)=Z(KK) 

10  CONTINUE 

C  WRITE  THE  INITIAL  VALUES  OF  THE  SEARCH. 

WRITE(4,500) 

500  F0RMAT(1X, 'THE  INITIAL  PERMUTATION  IS ' //5X, ' ELEMENT  NUMBER', 5X, 

1' POSITION' ) 

DO  600  11=1, IJ 

WRITE ( 4 , 5 0 1 )  RERR( RER( II) ) , IMR( II) 

501  F0RMAT(T1 1 , I3,T26 , 13) 

600  CONTINUE 
WRITE(4,507)  NACH 

507  FORMAT( IX, 'THE  NUMBER  OF  OBJECTIVES  IS',2X,I6) 

WRITE(4,502) 

502  FORMAT( IX, 'THE  INITIAL  ACHIEVEMENT  VECTOR  IS ’ //5X, ' PRIORITY  LEV 
1EL’ ,3X, 'ACHIEVEMENT' ) 

DO  601  11=1, NPRIOR 
WRITE( 4 , 503)  II,ZBEST(II) 

503  FORMAT(9X,I2, 13X.F10.4) 

601  CONTINUE 

C  BEGIN  THE  SEARCH  BY  CALLING  SUBROUTINE  PERMUTE. 

CALL  PERMUTE(& 100,&200,&300) 

C  WHEN  THE  SEARCH  ROUTINE  IS  TERMINATED,  CONTROL  RETURNS  TO  THE  MAIN 
C  ROUTINE  WITH  A  VALUE  OF  KK1 ,  INDICATING  THE  DEGREE  OF  SUCCESS. 

C  KK1=1  IS  WHEN  A  COMPLETE  CYCLE  OF  PERMUTATIONS  HAS  BEEN  COMPLETED 
C  WITHOUT  AN  IMPROVEMENT. 

100  KK1  =  1 


GOTO  110 

C  KK1=2  IS  WHEN  THE  NUMBER  OF  CYCLES  EQUALS  THE  MAXIMUM  NUMBER  OF 
C  ITERATIONS,  MAXITE. 

200  KK1=2 

GOTO  110 

C  KK1=3  IS  WHEN  THE  ACHIEVEMENT  VECTOR  EVALUATES  TO  ZERO,  e.g.  ALL 
C  THE  OBJECTIVES  HAVE  BEEN  COMPLETELY  SATISFIED. 

300  KK1=3 

110  CALL  WRITE(KKl) 

STOP 

END 

SUBROUTINE  PERMUTE( * ,* ,* ) 

C  THIS  ROUTINE  CONTROLS  THE  EXCHANGE  OF  ADJACENT  ELEMENTS  IN  THE 
C  PERMUTATION  ARRAY. 

C  K  IS  THE  PERMUTATION  INDEX  DURING  EACH  PASS  THROUGH  THE  PERMUTATION. 

C  ARRAY.  IEV=1  IS  THE  INDICATOR  OF  SUCCESS  OF  EACH  TRIAL  POINT. 

C  L  INDICATES  THE  NUMBER  OF  PASSES  THROUGH  THE  PERMUTATION  ARRAY. 

COMMON/ COMMO 3/  RER, IMR,RERR 

COMMON/ COMM04/  K,M,AK 

COMMON/ COMMO 5/  IEV 

COMMON/ COMMO 6/  MAXITE 

COMMON/ COMM1 2/  ZBEST 

COMMON/COMM15/  L,KK 

COMMON/ COMM1 6/  NPRINT, NCOUNT 

COMMON/COMM17/  MM,NPOOL,IJ 

COMMON/COMM18/  LSEED.ZLOLIM 

DIMENSION  S( 200) , IMR( 200) ,ZBEST( 10) 

INTEGER  RER( 200) , S , P,RERR( 200) , PI ,RERRT( 200) 

KK=1 

IEV=1 

L=1 

S(L)=0 

C  THE  FIRST  EXCHANGE 
1  P=RER(KK) 

RER(KK)=RER(KK+1 ) 

RER(KK+1)=P 

CALL  EVAL(&10,&200) 

10  RK=KK+ 1 

IF(IEV.EQ.l)  S(L)=S(L)+1 
IF(KK.EQ.IJ)  GOTO  100 
IF( IEV.EQ. 1 )  GOTO  50 

C  IF  THE  PREVIOUS  EXCHANGE  RESULTED  IN  NO  IMPROVEMENT,  THE  PREVIOUS  BASE 
C  POINT  IS  RETURNED  TO  BEFORE  THE  NEXT  EXCHANGE  IS  EFFECTED. 

P=RER(KK-1) 

RER(KK-1)=RER(KK) 

RER( KK) =RER( KK+ 1 ) 

RER(KK+1)=P 

CALL  EVAL(5il0,&200) 

C  IF  THE  PREVIOUS  EXCHANGE  RESULTED  IN  A  NEW  BASE  POINT,  THE  NEXT 
C  EXCHANGE  IS  FROM  THE  NEW  BASE  POINT. 


50  P=RER(KK) 

RER( KK) =RER( KK+ 1 ) 

RER(KK+1 )=P 

CALL  EVAL(&10,&200) 

C  IF  THE  EXCHANGE  BETWEEN  THE  LAST  TWO  ELEMENTS  RESULTS  IN  NO  NEW 
C  BASE  POINT,  THE  PREVIOUS  BASE  POINT  IS  RETURNED  TO. 

100  IF(IEV.EQ.l)  GOTO  120 
P=RER(KK-1) 

RER(KK-1 )=RER(KK) 

RER(KK)=P 
120  L=L+1 

IF(L.GT.MAXITE)  GOTO  175 
IF(S(L).EQ.S(L-1))  GOTO  150 
KK=1 
GOTO  1 

C  IF  A  COMPLETE  CYCLE  OF  PERMUTATIONS  HAS  NOT  BEEN  SUCCESSFUL 
150  IF(L.NE. 2)  GOTO  152 

C  IF  THE  PERMUTATIONS  FROM  THE  INITIAL  BASE  POINT  DO  NOT  RESULT  IN 
C  ANY  IMPROVEMENTS ; 

WRITE( 4,501) 

501  FORMAT( IX, 'THERE  IS  NO  IMPROVEMENT  ON  THE  INITIAL  PERMUTATION') 

C  IF  THE  BEST  SOLUTION  IS  NOT  LESS  THAN  ZLOLIM,  A  NEW  START  IS 

C  INTIATED.  THIS  ROUTINE  IS  CALLED  ONLY  6  TIMES.  THE  SIXTH 
C  CALL  CAUSES  THE  PROGRAM  TO  TERMINATE. 

152  IF(ZBEST(1).GT. ZLOLIM)  CALL  NEWSTART(&160) 

RETURN  1 

160  IF(NPRINT.NE. 1)  GOTO  503 
WRITE(4,502) 

502  FORMAT('OA  NEW  SEARCH  PATTERN  IS  NOW  BEING  USED.'/'  THE  ELEMENTS  A 
IRE  IN  THE  FOLLOWING  POSITIONS; '/5X, 'ELEMENT  NUMBER' , 5X, ' POSITION  N 
2UMBER' ) 

DO  503  11=1, IJ 

WRITE(4,504)  RERR(RER( II) ) , IMR( II) 

504  F0RMAT(T11,I3,T26,I3) 

503  CONTINUE 
KK=1 
GOTO  1 

C  IF  THE  MAXIMUM  NUMBER  OF  PERMUTATIONS  HAVE  BEEN  PERFORMED. 

175  RETURN  2 

C  IF  THE  ACHIEVEMENT  VECTOR  HAS  ZERO  VALUE  FOR  ALL  LEVELS. 

200  RETURN  3 
END 

SUBROUTINE  NEWSTART(*) 

C  THIS  SUBROUTINE  EFFECTS  A  NEW  START,  KEEPING  THE  ELEMENTS 
C  IN  THE  SAME  TRANSCUCER  ARRAY  POSITIONS,  BUT  REARRANGING  THE 
C  PERMUTATION  ARRAY. 

DIMENSION  IMR( 200) , NN( 200) , IMRR( 200) 

INTEGER  RER( 200) ,RERT( 200) 

COMMON/ COMMO 3/  RER, IMR.RERR 
COMMON/ COMMO 4/  K,M,AK 
COMMON/ COMM 16/  NPRINT , NCOUNT 


COMMON/ COMM1 7/  MM,NPOOL,IJ 

COMMON/COMM18/  LSEED,ZLOLIM 

DATA  NCOUNT/O/ 

IF(NCOUNT.GT.4)  RETURN 

NCOUNT=NCOUNT+ 1 

DO  10  11=1, IJ 

NN(II)=0 

CONTINUE 

LL=1 

ITER=0 

ITER=ITER+1 

Y=FLOAT( IJ)*RAN(LSEED)+1 .0 
NN1=INT(Y) 

IF(LL.NE.l)  GOTO  21 
NN(LL)=NN1 
LL=LL+1 
GOTO  1 

IFQTER.GT.  10*IJ)  GOTO  92 

DO  22  11=1, LL 

IF(NN1.EQ.NN(II) )  GOTO  l 

CONTINUE 

NN(LL)=NN1 

LL=LL+1 

IF(LL.GT.IJ)  GOTO  94 
GOTO  1 
ITER=5*IJ 
GOTO  1 

DO  30  11=1, IJ 
IMRR(II)=IMR(NN(II) ) 

RERT( 1I)=RER(NN( II) ) 
CONTINUE 
DO  40  11=1, IJ 
RER( II) =RERT (II) 

IMR( II)=IMRR( II) 

CONTINUE 
RETURN  1 
END 


SUBROUTINE  EVAL(*,*) 

THIS  ROUTINE  DETERMINES  IF  THE  CURRENT  TRIAL  POINT  GIVES  A  BETTER 
RESPONSE  THAN  THE  CURRENT  BASE  POINT.  IF  IT  DOES,  THE  TRIAL  POINT 
BECOMES  THE  NEW  BASE  POINT. 

THIS  ROUTINE,  WITH  ZVALUE,  AND  SOME  OF  THE  VARIABLE  NAMES 
ARE  FROM  CHARLES  ALVORD'S  DOCTORAL  THESIS. 

COMMO N/C0MM04/  K,M,AK 

COMMON/COMM05/  IEV 

COMMON/ COMM 1 1/  NACH , NPRIOR , NOBJ 

COMMON/COMM 12/  ZBEST 

COMMON/COMM1 5/  L,KK 

COMMON/ COMM 16/  NPRINT , NCOUNT 

DIMENSION  ZBEST( 10) , Z( 10) 

CALL  ZVALUE(Z) 
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58  CONTINUE 

C  EVALUATION  BY  COMPARISON  TO  THE  BEST  BASE  POINT  THUS  FAR. 

21  DO  20  JJ=1 , NPRIOR 

IF(Z(JJ).GT. ZBEST( JJ) )  GOTO  50 
IF(Z(JJ)+.Q001.LT.ZBEST( JJ) )  GOTO  30 
20  CONTINUE 

C  THE  TRIAL  POINT  RESULTS  IN  AN  EQUAL  REPONSE  TO  THE  CURRENT  BASE  POINT. 
IEV=0 
GOTO  50 

C  THE  TRIAL  POINT  RESULTS  IN  A  BETTER  RESPONSE. 

30  DO  35  JJ=1, NPRIOR 

IF(Z(JJ).NE.0.0)  GOTO  40 
35  CONTINUE 

C  IF  ALL  THE  OBJECTIVES  HAVE  BEEN  SATISFIED. 

GOTO  60 

C  INTIT IALIZE  THE  NEW  BASE  POINT. 

40  DO  45  JJ=1, NPRIOR 
ZBEST(JJ)=Z(JJ) 

45  CONTINUE 
IEV=1 

IF(NPRINT.NE. 1 )  GOTO  600 
WRITE(4,500)  L,KK 

500  FORMAT ( ' OA  NEW  BASE  POINT  HAS  BEEN  FOUND.'/'  THE  SUCCESSFUL  PERMUT 
1ATION  OCCURRED  DURING  CYCLE  //' , 12,'  AND  AT  PERMUTATION  POSITION', I 
23/ IX, 'THE  ACHIEVEMENT  VECTOR  IS’ //5X, 'PRIORITY  LEVEL' , 3X, 'ACHIEVEM 
3ENT' ) 

DO  600  I 1=1, NPRIOR 
WRITE(4,501)  II,ZBEST(II) 

501  FORMAT(8X, I4,8X,F10.4) 

600  CONTINUE 

RETURN  1 
50  IEV=0 

RETURN  1 

C  CONTROL  RETURNS  TO  PERMUTE  INDICATING  THAT  ALL  THE  OBJECTIVES  HAVE 
C  BEEN  SATISFIED. 

60  RETURN  2 
END 

SUBROUTINE  ZVALUE(Z) 

C  THIS  ROUTINE  DETERMINES  THE  LEVEL  OF  ACHIEVEMENT  OF  THE  OBJECTIVES. 
COMMON/ C0MM02/  Y 
COMMON/ COMM 1 1/  NACH, NPRIOR, NOB J 
COMMON/ COMM 12/  ZBEST 

COMMON/ C0MM1 3/  SIGN, IROW, WEIGHT, RHS.MPRI 

DIMENSION  Y( 3000) ,ZBEST( 10) , Z( 10) , SIGN( 2000) , IROW( 2000) 

DIMENSION  WEIGHT( 2000), RHS( 2000), MPRI( 2000) 

CALL  YVALUE 
DO  10  KK=1, NPRIOR 
Z(KK)=0.0 
10  CONTINUE 

DO  20  KK= 1 , NACH 
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IF(SIGN(KK) . LT.O.O)  GOTO  40 
DEV=Y( IROW( KK) )-RHS( IROW(KK) ) 

GOTO  50 

40  DEV=RHS( IROW(KK) )-Y( IROW(KK) ) 

50  IF(DEV. LT.O.O)  DEV=0.0 

Z(MPRI(KK) )=Z(MPRI(KK) )+WEIGHT(KK)*DEV 
20  CONTINUE 
RETURN 
END 

SUBROUTINE  DATAINN 

C  THIS  ROUTINE  READS  IN  THE  TRANSDUCER  ARRAY  PARAMETERS. 

C  THERE  ARE  2  MAIN  LOOPS  TO  READ  IN  THE  DATA,  ONE  FOR  SYMMETRIC 
C  AND  ONE  FOR  NON-SYMMETRIC 
C  IF  THE  ARRAY  IS  FOUR-FOLD  SYMMETRIC,  THERE 

C  ARE  2  ALTERNATE  WAYS  TO  INPUT  THE  DATA,  DEPENDING  ON  THE  EXISTENCE 
C  OF  AN  ARRAY  POSITION  AT  X=0,Y=0. 

COMMON/COMMOl/  IM , NEVEN , NROUTE 
COMMON/ COMM03/  RER, IMR,RERR 
COMMON/ C0MM04/  K,M,AK 
COMMON/ COMMO 6/  MAXITE 
COMMON/ COMMO 8/  C,D 
COMMON/ COMMIO/  A, AMP, ALPHA 
COMMON/ COMM1 4/  DOLAM.WAVE 
COMMON/COMM16 /  NPRINT , NCOUNT 
COMMON/ COMM1 7/  MM,NPOOL,IJ 
COMMON/ COMM18/  LSEED .ZLOLIM 

DIMENSION  A( 200) ,AMP( 200,200) ,ALPHA( 200, 200) ,C( 200) , D( 200) , IMR( 200) 
INTEGER  RER(200),RERR(200) 

DATA  D2RAD/ 1.74532925  E-02/ 

RAD2D=1 .0000/D2RAD 

READ(3,499)  DOLAM.WAVE 

READ( 3,500)  M .MAXITE , NROUTE , NPRINT , NPOOL 

IF(NPOOL.EQ. 1 )  READ( 3 , 506)  MM 

PI=4 . 0*DATAN( 1 . 0D00) 

AK=2*PI*DOLAM 

C  IF  NROUTE. NE. 1,  THE  ARRAY  IS  ASSYMMETRIC;  IF  NROUTE. EQ.l,  THE  ARRAY 
C  IS  SYMMETRIC. 

IJ=M 

IF ( NPOOL. EQ. 1 )  IJ=MM 
IK=*IJ 

IF( NROUTE. EQ. 1)  THEN 
READ(3,504)  NEVEN, IM 

C  IF  THE  ARRAY  IS  'ODD',  (i.e.  THERE  IS  AN  ELEMENT  AT  THE  CENTER  OF  THE 
C  ARRAY  X=0,Y=0) ,  THE  PARAMETERS  ARE  READ  IN  DIFFERENTLY  FOR  A  EVEN 
C  ARRAY.  NEVEN. EQ.l  IS  AN  ODD  ARRAY. 

IF(NEVEN.EQ. 1 )  THEN 
IK=IM+1 

C  COORDINATES  OF  THE  CENTER  ELEMENT  AND  THE  SHADING  VALUE  ARE  READ  IN 
C  FIRST. 

READ( 3 , 502)  C( 1) ,D( 1 ) ,A( 1) 

C  THE  COORDINATES  AND  SHADING  VALUES  FOR  THE  SYMMETRIC  POSITIONS  ARE 


C  READ  IN 


• 

DO  13  11=2, IK 
READ( 3, 502)  C(II) ,D(II) ,A(II) 

C  THE  ENTIRE  ARRAY  VALUES  ARE  INITIALIZED  USING  THE  EVEN  SYMMETRY. 

C( II+IM)=-D( II) 

D(II+IM)=C(II) 

C(II+2*IM)=-C(II) 

D(II+2*IM)=-D(II) 

C( II+3*IM)=D( II) 

D( II+3*IM)=-C( II) 

A(II+IM)=A(II) 

A( II+2*IM)=A( II) 

A( I 1+3* IM) =A( I I ) 

13  CONTINUE 

C  IF  THE  ARRAY  IS  EVEN. 

C  THE  COORDINATES  AND  SHADING  VALUES  OF  THE  SYMMETRIC  POSITIONS  ARE 

C  READ  IN.  THEN  THE  ENTIRE  ARRAY  VALUES  ARE  INITIALIZED  ACCORDING 
C  TO  THE  EVEN  SYMMETRY. 

ELSE 

1 2  IK=IM 

DO  17  11=1, IK 

READ(3,502)  C(II) ,D(II) ,A(II) 

C( II+IM)=-C( II) 

C( II+2*IM)=-C( II) 

C( II+3*IM)=C( II) 

D( II+IM)=D( II) 

D( II+2*IM)=-D( II) 

D(II+3*IM)=-D( II) 

A(  II+IM)=A( II) 

A( I 1+2* IM) =A( II) 

A(II+3*IM)=A(II) 

17  CONTINUE 

ENDIF 

C  I?  THE  ARRAY  IS  NOT  4-FOLD  SYMMETRIC  THE  COORDINATES,  AND  SHADING 
C  VALUES  OF  THE  ENTIRE  ARRAY  ARE  READ  IN. 


11 

ELSE 

DO  10  11=1, 

IJ 

READ(3,502) 

C(II),D(II),A(II) 

10 

CONTINUE 

ENDIF 

C  IF  THERE  ARE  ELEMENTS  IN  THE  'POOL’,  THEIR  EQUIVALENT  ARRAY  POSITION 
C  AND  SHADING  VALUES  MUST  ALSO  BE  SPECIFIED,  (i.e.  SET  EQUAL  TO  ZERO.) 
14  IF ( NPOOL . NE . 1 )  GOTO  22 

DO  21  II=M+1 ,MM 
C(II)=0.0 
D(II)=0.0 
A(II)=0.0 
21  CONTINUE 

C  ELEMENT  TOLERANCES  ARE  READ  IN  IDENTICALLY  FOR  ALL  LOOPS. 

C  THE  OUTER  LOOP  CONTROLS  INPUT  ACCORDING  TO  ELEMENT,  THE  INNER 
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C  ACCORDING  TO  POSITION. 

22  IF(NROUTE.NE. 1 )  THEN 
DO  20  11=1, IJ 
READ(3,506)  RERR(II) 

DO  20  JJ=1 , IK 

READ( 3 , 503)  AMP( II , JJ) , ALPHA 1 
C  CONVERT  TOLERANCES  FROM  DEGREES  TO  RADIANS. 

ALPHA( I I , J J ) =ALPHA 1 *D2 RAD 
20  CONTINUE 
ELSE 

DO  27  11=1, IJ 
READ( 3,506)  RERR(II) 

DO  28  JJ=1 , IK 

READ(3,503)  AMP(II,JJ) .ALPHA  1 
ALPHAC I I , J J ) =ALPHA 1 * D2RAD 

C  SINCE  THE  RESPONSE  OF  AN  ELEMENT  WILL  BE  THE  SAME  FOR 
C  ANY  OF  THE  FOUR  SYMMETRIC  POSITIONS,  ARRAY  POSITION  MUST  BE 
C  INITIALIZED. 

IF(NEVEN.EQ. l.AND. JJ.EQ. 1)  GOTO  28 
AMP( II , JJ+IM)=AMP( II, JJ) 

AMP( II , JJ+2*IM)=AMP( II , JJ) 

AMP( I I , J J+3* IM) =AMP( I I , J J ) 

ALPHA( I I , J J+IM) =ALPHA( I I , J J ) 

ALPHA( I I , J J+2* IM) =ALPHA( I I , J J ) 

ALPHAC I I , J J+3* IM) = ALPHA( I I , J J ) 

28  CONTINUE 

DO  27  KK=4*IM+1 , MM 
AMPC I I , KK) =0 . 0 
ALPHAC I I, KK) =0.0 
27  CONTINUE 

ENDIF 

C  READ  IN  THE  STARTING  POSITIONS  FOR  EACH  ELEMENT. 

DO  40  JJ=1,IJ 

READC  3 , 504)  IMRC JJ) ,RERC JJ) 

40  CONTINUE 
42  CALL  GOALIN 

READC 3, 505)  LSEED.ZLOLIM 
CALL  NEARIN 

499  FORMAT( 2F9. 5) 

500  FORMATC5I5) 

502  FORMATC3F9.5) 

503  FORMAT(2F9.5) 

504  FORMATC  215) 

505  FORMAT(I10,F9.4) 

506  FORMATC 15) 

RETURN 

END 

SUBROUTINE  GOALIN 

C  THIS  ROUTINE  READS  IN  MATH  PROGRAMMING  PARAMETERS,  SUCH  AS  THE 
C  THETA  AND  PHI  ANGLES  THAT  ARE  CONSTRAINED. 
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C  THE  GOAL  PROGRAMMING  PARAMETERS,  SUCH  AS  THE  RHS  VALUES,  THE 
C  WEIGHTS,  AND  SIGN  OF  THE  DEVIATION  VARIABLES  ARE  ALSO  INPUT 
C  HERE. 

COMMON/ COMMO 7/  THETA, DELTHE 

COMMON/ C0MM09/  PHI,NPHI 

COMMON/COMM1 1/  NACH,NPRIOR,NOBJ 

COMMON/COMM13/  SIGN, IROW, WEIGHT, RHS, MPRI 

DIMENSION  RHS( 2000) ,SIGN( 2000) , IROW( 2000) , WEIGHT/ 2000) 

DIMENSION  MPRI( 2000) ,THETA( 100) ,PHI(40) 

DATA  RAD2D/57. 2958/, D2RAD/1. 74532925  E-02/ 

C  READ  IN  THE  PHI  ANGLES  CONSTRAINED. 

READ(3,99)  NREAD1 
99  FORMAT/ 13) 

IF/NREAD1.EQ. 1)  THEN 
READ/3,99)  NPHI 
DO  103  11=1, NPHI 
READ/3,101)  PHI2 
PHI/ II)=PHI2*D2RAD 
103  CONTINUE 
ELSE 

READ/ 3, 100)  PHI 1, DELPHI 
TEM=Z  360. 0-PHI1) /DELPHI 
NPHI=TEM 
TEMl=INT/TEM)/2 

IF/TEM/2.0.NE.TEM1)  NPHI=TEM+1.0 
PHI/ 1)=PHI1*D2RAD 
DELPHI=DF  ,PHI*D2RAD 
DO  11  11=2, NPHI 
PHI/ II)=PHI/ 1 )+( II-l )*DELPHI 
1 1  CONTINUE 
ENDIF 

READ/ 3, 101)  DELTHE 
NACH=0 

C  READ  IN  THE  THETA  VALUES  CONSTRAINED  FOR  EACH  PHI. 

DO  10  11=1, NPHI 
READ/ 3, 101)  THETA/ II) 

C  SUM  OVER  ALL  ANGLES  CONSTRAINED  TO  DETERMINE  THE  NUMBER 
C  OF  OBJECTIVES. 

NACH=NACH+/ 90. O-THETA/ II) )/DELTHE+l 
THETA/ 1 1 ) =THETA/ 1 1 ) *  D2  RAD 
10  CONTINUE 

C  CONVERT  DEGREES  TO  RADIANS. 

DELTHE=DELTHE*D2RAD 

NOBJ=NACH 

READ/ 3, 102)  NSTED 

C  /NSTED. EQ.l)  MEANS  THAT  SIGN, WEIGHT, MPRI , RHS  ARE  ALL  THE  SAME  FOR  EACH 
C  OBJECTIVE. 

IF/ NSTED. NE. 1)  THEN 

C  READ  IN  THE  VALUES  FOR  EACH  OBJECTIVE. 

READ/ 3, 102)  NPRIOR 
DO  30  I 1=1, NOB J 

READ/ 3, 104)  MPRI/ 1 1), WEIGHT/ 1 1) , IROW/ II), SIGN/ 1 1), RHS/ II) 
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)  CONTINUE 
ELSE 

READ  IN  THE  SIDE  LOBE  LEVEL  AND  SPECIFY  ALL  THE  OTHER  VALUES 
NPRIOR= 1 

READ( 3,101)  RHS 1 
DO  60  JJ=1 , NOBJ 
MPRI( J J)= 1 
WEIGHT( J J)= 1 
SIGN( JJ)=-1 . 0 
RHS( JJ)=RHS1 
IROW( JJ)=JJ 
)  CONTINUE 
ENDIF 

DO  FORMAT( 2F8.4) 

D1  FORMAT(F8.4) 

D2  FORMAT (12) 

D4  FORMAT(I2,F8.4,I4,2F8.4) 

RETURN 

END 

SUBROUTINE  WRITE(KKI) 

THIS  ROUTINE  IS  CALLED  ONLY  AFTER  THE  SEARCH  HAS  BEEN  TERMINATED.  IT 
WRITES  OUT  ALL  PERTINENT  DATA  BASED  ON  THE  RESULTS  OF  THE  SEARCH. 
COMMON/ C0MM01/  IM,NEVEN,NROUTE 
COMMON/ C0MM03/  R£R,IMR,RERR 
COMMON/ C0MM04/  K , M , AK 
COMMON/ COMM05/  IEV 
COMMON/ COMM 10/  A, AMP, ALPHA 
COMMON/ C0MM11/  NACH,NPRIOR,NOBJ 
COMMON/ C0MM1 2/  ZBEST 
COMMON/ C0MM15/  L,KK 
COMMON/ C0MM1 6/  NPRINT , NCOUNT 
C0MM0N/C0MM17/  MM.NPOOL.IJ 

DIMENSION  IMR( 200) ,AMP( 200, 200) ,ALPHA( 200, 200) ,ZBEST( 10) ,A(200) 
INTEGER  RER( 200) , RERR( 200) 

DATA  D2RAD/ 1.74532925  E-02/ 

RAD2D=l.0000/D2RAD 

KK1  INDICATES  THE  SUCCESS  OF  THE  SEARCH. 

THE  RESULTS  ARE  PRINTED  ACCORDING  TO  THE  SUCCESS  OF  THE  SEARCH. 

GOTO( 10,20, 30) ,  KK1 

IF  A  COMPLETE  CYCLE  OF  PERMUTATIONS  HAS  BEEN  COMPLETED  AND 
NO  IMPROVEMENTS  WERE  MADE. 

0  WRITE( 4 , 100) 

DO  FORMAT (IX, 'A  COMPLETE  CYCLE  OF  PERMUTATIONS  HAS  BEEN  MADE  WITHOUT 
1  ANY' /IX, 'IMPROVEMENTS  ON  THE  PREVIOUS  BASE  POINT.'/) 

GOTO  50 

0  WRITE( 4 , 200) 

00  FORMAT( IX, 'THE  MAXIMUM  NUMBER  OF  CYCLES  OF  THE  PERMUTATION  ALGORITH 
1M  HAVE  BEEN  COMPLETED.') 

GOTO  50 

0  WRITE( 4 , 300) 

00  FORMAT ( ’  THE  ACHIEVEMENT  FUNCTION  IS  0.0.  ALL  THE  OBJECTIVES  HAVE 
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1BEEN  SATISFIED'/) 

GOTO  60 

50  IF ( NPRI NT . EQ . 1 )  GOTO  51 
WRITE(4, ill)  L , NCOUNT 

111  FORHATC  THE  NUMBER  OF  CYCLES  COMPLETED  IS  ',12,'.  THE  NUMBER  OF 
1NEW  SEARCH  PATTERNS  USED  IS  ',12,'.') 

51  WRITE(4 , 101 ) 

101  FORMAT( IX, 'THE  ACHIEVEMENT  FUNCTION  IS ' /T1 5 PRIORITY  LEVEL', 5X,' A 
1CHIEVEMENT' ) 

DO  1  11=1 , NPRIOR 

WRITE( 4 , 102)  II.ZBEST(II) 

102  FORMAT(T17 ,I2,T33,F12.2) 

1  CONTINUE 

60  WRITE(4 , 103) 

103  F0RMAT(1X, 'THE  ELEMENTS  ARE  IN  THE  FOLLOWING  POSITIONS ' /5X, ' ELEMEN 
IT  NUMBER' ,5X, 'POSITION  NUMBER') 

DO  2  11=1, IJ 

WRITE (  4 , 104)  RERR(  RER(  II)  ) ,  IMR(  II) 

104  FORMAT(T 1 1,I3,T30,I3) 

2  CONTINUE 

r* 

C  THIS  WRITES  THE  ELEMENT  TOLERANCES  INTO  THE  FILE  CREATED  BY  THE 
C  SUBROUTINE  NEARIN,  SO  THAT  NEARPRES  CAN  READ  THAT  FILE  DURING 
C  EXECUTION.  THE  PORTION  OF  THE  FILE  WRITTEN  HERE  IS  THE  BEST  SOLUTION 
C  FOUND  BY  THE  PERMUTATION  SEARCH  ALGORITHM. 

C 

500  DO  601  11=1, IJ 

WRITE( 2 , 602)  IMR(II) 

601  CONTINUE 

602  F0RMAT(I4) 

DO  600  11=1, IJ 

WRITE( 2,700)  AMP(RER( II) , IMR( II) ) , ALPHA( RER( II) , IMR( II) )*RAD2D 
700  F0RMAT(F9. 5 , ' , ' ,F9. 5) 

600  CONTINUE 
900  WRITE( 2,403) 

403  FORMAT('O') 

RETURN 

END 

SUBROUTINE  YVALUE 

C  THIS  ROUTINE  DETERMINES  THE  RESPONSE  OF  THE  ARRAY  FOR  EACH  THETA  AND 
C  PHI  THAT  IS  CONSTRAINED,  AS  MEASURED  AGAINST  THE  MAIN  LOBE  RESONSE. 
COMMON/ COMMO 2/  Y 
COMMON/ COMMO 3/  RER, IMR,RERR 
COMMON/ C0MM04/  K,M,AK 
COMMON/ COMMO 7/  THETA, DELTHE 
COMMON/ COMM08/  C,D 
COMMON/ COMM09/  PHI,NPHI 
COMMON/COMMlO/  A, AMP, ALPHA 
COMMON/COMM17/  MM,NPOOL,IJ 
COMMON/COMM19/  DNORM 
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DIMENSION  IMR( 200) ,Y( 3000) ,C( 200) , D( 200) ,A(200) ,AMP( 200 , 200) 
DIMENSION  ALPHA( 200, 200) ,PHI(40) ,THETA( 100) 

INTEGER  RER( 200) 

DATA  RAD2D/57. 2958/, D2RAD/1. 74532925  E-02/ 

REAL* 8  SUMR , SUMI , DCOS , DSIN , U 

C  ZZ1  REPRESENTS  THE  RESPONSE  OF  THE  MAIN  BEAM.  ZZ  IS  THE  RESPONSE  AT  A 
C  SPECIFIC  THETA  AND  PHI.  ZZ1-ZZ  REPRESENTS  THE  DIFFERENCE  (IN  dB  ) 

C  BETWEEN  THEM. 

CALL  INITIALIZE(ZZl) 

KK=0 

1=1 

DO  50  JJ=1 , NPHI 

C  EACH  THETA(I)  REPRESENTS  THE  THETA  ANGLE  WHERE  THE  SIDE  LOBE  STARTS. 

C  FOR  EACH1  PHI  CONSTRAINED. 

THETAl=THETA( I) 

11  SUMR=0.0 
SUMI=0.0 

USIN=SIN(THETA1 )*SIN( PHI( JJ) ) 

UCO S=SIN( THETA 1 )*COS( PHI( JJ) ) 

C  SUM  OVER  ALL  THE  ELEMENTS  TO  DETERMINE  THE  RESPONSE  OF  THE  ARRAY. 

DO  20  11=1, IJ 

IF ( A( IMR( II)).EQ.0.0)  GOTO  20 
U1=C( IMR( II) )*UC0S 
U2=D( IMR( II) )*USIN 

U=DBLE( AK*(U1+U2 )+ALPHA(RER( II) , IMR( II) ) ) 

H=A( IMR( II) )*AMP(RER( II) , IMR( II) ) 

SUMR=SUMR+H*DCOS(U) 

SUMI=SUMI+H*DSIN(U) 

20  CONTINUE 
KK=KK+1 

AMPL=SNGL( SUMR* SUMR+SUMI* SUMI) 

AMPL=SQRT(AMPL) 

SAMP=AMPL/ DNORM 

IF( ABS( SAMP) . GT. 1 .0592E-5)  THEN 
ZZ=20*ALOGI0(SAMP) 

ELSE 

13  ZZ=-99.0 
ENDIF 

C  SUBTRACT  THE  SIDE  LOBE  RESPONSE  FROM  THE  MAIN  BEAM  RESPONSE. 

C  THIS  GIVES  THE  VALUE  OF  THE  LEFT  HAND  SIDE  OF  THE  OBJECTIVE. 

14  Y(KK)=ZZ1-ZZ 

C  INCREMENT  THE  THETA  ANGLE. 

THETA 1=THETA1+DELTHE 

C  IF  THETA  EXCEEDS  90.0  DEGREES,  INCREMENT  THE  PHI  ANGLE. 
IF(THETA1*RAD2D.GT.90.0)  GOTO  45 
GOTO  11 

C  INCREMENT  THE  PHI  ANGLE  AND  THE  THETA  ANGLE  INDEX,  I. 

45  I-I+l 
50  CONTINUE 

C  IF  THE  PHI  ANGLE  EXCEEDS  360.0  DEGREES,  THE  RESPONSE  AT  ALL 
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CONSTRAINED  ANGLES  HAS  BEEN  DETERMINED. 
RETURN 
END 


SUBROUTINE  INITIALIZE(ZZl) 

THIS  SUBROUTINE  DETERMINES  THE  RESPONSE  OF  THE  ARRAY  ON  THE  MAINBEAM 
i.e.  THETA  =0.0,  AND  PHI  =0.0 
COMMON/ C0MM03/  RER, IMR,RERR 
COMMON/ COMM04/  K,M,AK 
COMMON/ C0MM0 8/  C,D 
COMMON/ COMM 10/  A, AMP, ALPHA 
COMMON/ C0MM17/  MM,NPOOL,IJ 
COMMON/ COMM1 9/  DNORM 

DIMENSION  IMR( 200) ,C( 200) , D( 200) ,A( 200) ,AMP( 200,200) 

DIMENSION  ALPHA(200,200) 

INTEGER  RER( 200) 

DATA  RAD2D/57. 2958/, D2RAD/1. 74532925  E-02/ 

DATA  INDEX/0/ 

REAL* 8  SUMR, SUMI , DCOS , DSIN ,U 
KK=  1 

IF ( INDEX. GT. 0 )  GOTO  23 
DNORM=0. 0 
DO  23  11=1, M 
DNORM=DNORM+A( II) 

INDEX=INDEX+1 
3  CONTINUE 
THETA2=0. 0 
PHI2=0.0 
SUMR=0 . 0 
SUMI=0.0 

USIN=SIN(THETA2)*SIN( PH 12) 

UCOS=SIN(THETA2)*COS( PHI2) 

SUM  OVER  ALL  THE  ELEMENTS  TO  DETERMINE  THE  RESPONSE  OF  THE  ARRAY. 

DO  20  11=1, IJ 

IF( A( IMR( II) ) . EQ. 0. 0)  GOTO  20 
U1=C( IMR( II) )*UCOS 
U2=D( IMR( II))*USIN 

U=DBLE ( AK* ( U 1+U2 ) +ALPHA( RER( II) , IMR( II) ) ) 

H=A( IMR( II) )*AMP( RER( II) , IMR( II)) 

SUMR=SUMR+H*DCOS(U) 

SUMI=SUMI+H*DSIN(U) 

3  CONTINUE 

AMPL=SNGL( SUMR*SUMR+SUMI*SUMI) 

AMPL=SQRT(AMPL) 

SAMP=AMPL/DNORM 

IF( ABS( SAMP) .GT. 1.0592E-5)  THEN 
ZZ1=20*ALOG10(SAMP) 

ELSE 

3  ZZl=-99.0 

ENDIF 
RETURN 
END 
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SUBROUTINE  NEARIN 

C  NEARIN  WRITES  THE  TRANSCUCER  ARRAY  PARAMETERS  IN  AN  OUTPUT  FILE  IN  A 
C  FORMAT  WHICH  IS  COMPATIBLE  WITH  THE  INPUT  FORMAT  FOR  THE  NEARPRES 
C  PROGRAM.  THIS  ROUTINE  CONTROLS  ONLY  THE  INITIAL  PARAMETERS. 

C  THE  PARAMETERS  OF  THE  FINAL  SOLUTION  ARE  PRINTED  BY  THE  WRITE 
C  SUBROUTINE  THIS  ROUTINE  INCORPORATES  THE  SYMMETRY  OF  THE  ARRAY  INTO 
C  THE  OUTPUT  FILE. 

COMMON/COMMOl/  IM , NEVEN , NROUTE 
COMMON/ COMM03/  RER, IMR,RERR 
COMMON/COMM04/  K,M,AK 
COMMON/ COMM08/  C,D 
COMMON/ COMM 10/  A, AMP, ALPHA 
COMMON/ C0MM1 4/  DOLAM.WAVE 
C0MM0N/C0MM17/  MM,NPOOL,IJ 

DIMENSION  IMR( 200) ,A(200) ,AMP( 200,200) ,ALPHA( 200,200) ,C( 200) ,D( 200) 
INTEGER  RER( 200) 

DATA  D2RAD/ 1.74532925  E-02/ 

RAD2D=1 . 0000/D2RAD 
IF(NPOOL.EQ.l)  THEN 
WRITE( 2,400)  NPOOL.MM 

400  FORMATC  1,1,1' / '-90.0,90.0,2.0'/ '0.0, 180.0,5.0  712/13) 

ELSE 

WRITE(2,401)  NPOOL 

401  FORMATC 1,1,1' /' -90. 0,90. 0,2. O' /' 0.0, 180. 0,5. O’ /I2) 

ENDIF 

500  WRITE( 2,711)  NROUTE, M 

711  FORMATC  2,  2 , '  ,  12 , '  , '  ,  12 , '  ,  1') 

C  WRITE  THE  ARRAY  PARAMETERS  USING  FOUR-FOLD  SYMMETRY. 

IF( NROUTE. EQ. 1)  THEN 

N1L=INT( SQRT(FLOAT(IM) ) ) 

IF(NEVEN.EQ.l)  THEN 
N2L=N1L+1 
ELSE 

714  N2L=N1L 

ENDIF 

719  WRITE(2,712)  IM, NEVEN, NIL, N2L 

712  FORMATC 12,' , ' ,12, ' , ' ,12, ' , ' ,12) 

IF( NEVEN. ME. 1)  THEN 

C  WRITE  THE  ARRAY  PARAMETERS  USING  EVEN  SYMMETRY. 

C  WRITE  THE  SHADING 

DO  601  11=1, IM 
WRITE( 2 ,702)  A(II) 

702  FORMATC F8. 4) 

601  CONTINUE 

C  WRITE  THE  ELEMENT  TOLERANCES. 

901  DO  902  11=1, IJ 
WRITEC  2 , 302)  IMRCII) 

902  CONTINUE 

DO  603  11=1, IJ 

WRITEC2.700)  AMPCRER(II) , IMR( II) ) ,ALPHA(RER( II) , IMR( II) )*RAD 


12D 


603  CONTINUE 

700  FORMAT(F9.5, ' , ' ,F9.5) 

C  WRITE  THE  COORDINATES. 

DO  599  11=1, IM 
WRITE(2,701)  C( II) , D( II) 

701  F0RMAT(F6. 3 , ' , ' ,F6. 3) 

599  CONTINUE 

WRITE( 2,725)  WAVE , DOLAM 
725  F0RMAT(F8. 4/F8. 4) 

C  WRITE  THE  ARRAY  PARAMETERS  USING  ODD  SYMMETRY. 

ELSE 

710  DO  721  11=1, IM+1 

WRITE(2,702)  A(II) 

721  CONTINUE 

DO  903  11=1, IJ 
WRITE(2,302)  IMR(II) 

903  CONTINUE 

302  FORMAT( 14) 

DO  724  11=1, IJ 

WRITE(2,700)  AMP( RER( II) , IMR( II) ) , ALPHA(RER( II) , IMR( II) )*RAD 

12D 

724  CONTINUE 

C  WRITE  THE  COORDINATES. 

DO  723  11=1, IM+1 
WRITE( 2,701)  C(II) ,D(II) 

723  CONTINUE 

WRITE( 2 ,725)  WAVE, DOLAM 
ENDIF 
ELSE 

C  WRITE  THE  PARAMETERS  FOR  A  NON-SYMMETRIC  ARRAY. 

DO  801  II=l,IJ 
WRITE(2,702)  A(II) 

801  CONTINUE 

DO  909  11=1, IJ 
WRITE(2,302)  IMR(II) 

909  CONTINUE 

DO  904  11=1, IJ 

WRITE(2,700)  AMP(RER( II) , IMR( II) ) , ALPHA( RER( II) , IMR( II) )*RAD 

12D 

904  CONTINUE 

DO  803  11=1, IJ 
WRITE(2,701)  C( I I) , D( I I) 

803  CONTINUE 

WRITE( 2,804)  WAVE 

804  FORMAT(F8.4) 

ENDIF 

900  RETURN 
END 
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APPENDIX  E 

LISTING  OF  NEO-KENDIG  COMPUTER  SELECTION  PROGRAM 


C********************************************************************** 

H  K  V 

COMPUTER  IMPLEMENTATION  OF  THE  NEO-KENDIG  SELECTION  TECHNIQUE 
WRITTEN  BY  KING  W.  WIEMANN 

C************************************* **************** ***************** 

INTEGER  U( 200,2) 

REAL  A( 2) , V( 2) 

COMMON/ COMMOl/  M,MM,IM,IMl 

COMMON/COMM02/  U 

COMMON/ C0MM07/  A,V 

COMMON/ C0MM1 3/  NEVEN,N1L,N2L,LSEED 

CALL  DATAININ 

CALL  RESPONSE(l) 

CALL  STAT(l) 

CALL  DELETE 
CALL  STAT(l) 

WRITE( 4 , 500)  100.*A( 1 ) ,A( 2) 

500  F0RMAT(5X, 'AMPLITUDE  MEAN  =  ',F12.8,'  PHASE  MEAN  =  ',F12. 
1) 

WRITE(4,501)  100.*V( 1 ) ,V( 2) 

501  F0RMAT(5X, 'AMPLITUDE  VARIANCE  =  ’,F12.8,’  PHASE  VARIANCE  =  ’  ,F12, 

l/) 

WRITE(4,502) 

502  FORMAT(90X, 'GROUP' ,5X, ' POSITION  : ' ,3X, ’ 1 ’ ,5X, ’ 2' ,5X, ' 3 ' ,5X, ' 4 ' ) 
DO  10  11=1, IM1 

IF  (II.GT.l)  CALL  RESPONSE(II) 

CALL  REMOVE ( II) 

IF(II.EQ. l.AND.NEVEN.EQ. 1)  GOTO  10 
CALL  OPTIMO(II) 

CALL  WRITEIT(II) 

10  CONTINUE 

CALL  WRITEIT(O) 

CALL  NEARIN 1 
CALL  PERMSERIN 
STOP 
END 

SUBROUTINE  CMAX(D1,S) 

C  ROUTINE  TO  DETERMINE  THE  MAXIMUM  VALUE  OF  C.  USED  FOR  SORTING  THE 
C  REMAINING  ELEMENTS. 

INTEGER  S 
DIMENSION  C( 20) 

COMMON/ C0MM17/  C 
Dl=l . OE-1 2 
DO  10  LL= 1 , S 

IF  (C(LL).GT.Dl)  D1=C(LL) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  DATAININ 

C  ROUTINE  TO  READ  IN  THE  PARAMETERS  NECESSARY  TO  IMPLEMENT  THE 
C  NEO-KENDIG  TECHNIQUE. 

REAL* I 6  S(200,9) 

INTEGER  U(200,2),OSS(40) ,TERM(3) ,RER(200) ,RERR(200) , IMR(200) 
REAL  RA(40) ,XA(40) ,X(70) , Y(70) ,SHADE(70) 

CHARACTER* 1  NAME (40) 

COMMON/COMMOl/  M,MM,IM,1M1 

COMMON/ COMMO 2/  U 

COMMON/ C0MM06/  OSS, TERM 

COMMON/ C0MM6A/  NAME 

COMMON/ COMMO 8/  RA,XA,F,SIZE 

COMMON/ C0MM8 A/  S 

COMMON/ C0MM12/  RER,RERR, IMR 

COMMON/ C0MM1 3/  NEVEN,N1L,N2L,LSEED 

COMMON/ COMM1 4/  WAVE , DOLAM 

COMMON/ COMM1 5/  X,Y, SHADE 

READ(3,500)  MM,M,IM,IM1 ,NEVEN 

500  FORMAT( 515) 

C  OSS  IS  THE  ORDER  OF  SEARCH.  IT  SPECIFIES  IN  WHICH  ORDER  THE  ARRAY 
C  IS  TO  BE  FILLED. 

READ(3,503)  (OSS( JJ) , JJ=1 ,IM1) 

503  FORMAT( 2013) 

C  TERM  CONTROLS  WHICH  KENDIG  TERMS  ARE  TO  BE  THE  CRITERIA. 

C  TERMS  3,  4  AND  8  ARE  ALLOWED  IN  ANY  GROUPING  OF  1,  2,  OR  3  TERMS. 
READ(3,505)  (TERM( JJ) , JJ=1 ,3) 

505  FORMAT( 312) 

READ( 3,501)  F.SIZE 

501  FORMAT(2F10.4) 

READ( 3 , 508)  LSEED 

508  FORMAT( 19) 

C  U  IS  THE  NUMBER  ASSIGNED  TO  EACH  ELEMENT;  S  CONTAINS  THE 
C  EQUIVALENT  CIRCUIT  VALUES. 

DO  10  JJ=1 ,MM 
U( JJ,2)=0 

READ(3,502)  U( JJ, 1 ) , (S( JJ,LL) ,LL=2 , 7 ) 

502  FORMATQ5.6D15.9) 

10  CONTINUE 

C  RA,  AND  XA  ARE  THE  REAL  AND  IMAGINARY  RADIATION  LOADINGS  FOR 
C  EACH  POSITION.  NAME  IS  THE  LETTER  ASSIGNED  TO  EACH  SYMMETRIC 
C  POSITION. 

DO  30  JJ=1 , IM1 

READ( 3 , 504)  RA(OSS( JJ) ) , XA(OSS( JJ) ) , NAME(OSS( JJ) ) 

504  F0RMAT(2E12.4,A1) 

30  CONTINUE 

C  NIL  AND  N2L  SPECIFY  THE  LENGTH  AND  WIDTH  OF  THE  FIRST  QUADRANT, 

C  (IN  NUMBER  OF  ELEMENTS). 

NIL- INT(SQRT(FLOAT( IM) ) ) 

IF(NEVEN.EQ. 1 )  THEN 
N2L=NIL+1 
ELSE 
N2L=N1L 
ENDIF 


C  RER  IS  THE  SEQUENTIAL  NUMBER  OF  EACH  ELEMENT  OF  THE  ARRAY  IN  ARRAY 
C  ORDER. 

C  RERR  IS  THE  ASSIGNED  NUMBER  OF  EACH  ELEMENT  OF  THE  ARRAY  IN  ARRAY 
C  ORDER. 

C  IMR  IS  THE  ORDER  OF  TRANSDUCER  ARRAY  POSITIONS  IN  THE  PERMUTATION 
C  ARRAY. 

DO  40  KK= 1 , MM 
RER(KK)=0 
RERR(KK)=0 
IMR(KK)=KK 
40  CONTINUE 
N=0 

IF(NEVEN.EQ. 1)  N=1 

C  X  AND  Y  ARE  THE  X  AND  Y  COORDINATES  OF  EACH  SYMMETRIC  POSITION, 

C  AND  SHADE  IS  THE  SHADING  VALUE  FOR  EACH  POSITION. 

DO  SO  JJ=I , IM+N 

50  READ( 3 , 507)  X( JJ) , Y( JJ) , SHADE( JJ) 

507  FORMAT ( 3F8. 4) 

READ(3,508)  WAVE , DOLAM 

508  FORMAT( 2F9. 4 ) 

RETURN 

END 

SUBROUTINE  DELETE 

C  ROUTINE  TO  PERMANENTLY  DELETE  ELEMENTS  FROM  THE  SELECTION  POOL  BEFORE 
C  THE  SELECTION  PROCESS  BEGINS. 

INTEGER  U(200,2) 

COMPLEX* 16  ELEMNT( 200) 

REAL  A( 2) , V( 2) 

COMMON/COMMOl/  M,MM,IM,IM1 
COMMON/ COMMO 2/  U 
COMMON/ COMMO 7/  A,V 
COMMON/ COMM1 6/  ELEMNT 
DO  10  JJ=1 ,MM 
T1=DREAL(ELEMNT( JJ) )-A( 1 ) 

T2=DIMAG(ELEMNT( JJ) )-A( 2) 

IF  (Tl. LT. 4. 0*V( 1 ) )  GOTO  10 
IF  (T2.LT.4.0*V(2))  GOTO  10 
WRITE( 4 , 100)  U(JJ,1) 

100  FORMAT ( *  ELEMENT  NUMBER  ',13,'  WAS  DELETED  FROM  CONSIDERATION') 
U(JJ,1)=0 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  DMAX(Dl.S) 

C  ROUTINE  TO  FIND  THE  MAXIMUM  VALUE  OF  THE  D  ARRAY. 

C  USED  TO  FIND  THE  BEST  PERMUTATIONS  THAT  MEET  THE  KENDIG  TERM 
C  UNDER  CONSIDERATION. 

INTEGER  S 
DIMENSION  D( 40) 

COMMON/ C0MM10/  D 
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D1=1.0E-12 
DO  10  KK=1,S 

IF  (D(KK).GT.Dl)  D1=D(KK) 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  NEARIN1 

C  NEARIN 1  WRITES  THE  TRANSCUCER  ARRAY  PARAMETERS  IN  AN  OUTPUT  FILE  IN  A 
C  FORMAT  WHICH  IS  COMPATIBLE  WITH  THE  INPUT  FORMAT  FOR  THE  NEARPRES 
C  PROGRAM. 

REAL  X(70),Y(70) ,SHADE(70) 

INTEGER  U(200,2) ,RER(200) 

COMPLEX* 16  ELEMNT( 200) 

COMMON/ COMMO l /  M,MM,IM,IMl 
COMMON/COMM02/  U 
COMMON/ C0MM1 2/  RER,RERR, IMR 
C0MM0N/C0MM13/  NEVEN, NIL, N2L.LSEED 
C0MM0N/C0MM14/  WAVE , DOLAM 
COMMON/ C0MM1 5/  X,Y, SHADE 
COMMON/ COMM 16/  ELEMNT 
WRITE(2,500)  MM 

500  FORMAT( ,0,l,l,/'-90.0,90.0,2.0,/,0.0,180.0,5.0,/,l,/I3) 

WRITE( 2,501)  M 

501  FORMATC 1,0,1,’ ,12,' ,1') 

C  WRITE  THE  ARRAY  PARAMETERS  USING  FOUR-FOLD  SYMMETRY. 

WRITE( 2 , 502)  IM , NEVEN , N 1 L , N2L 

502  F0RMAT(I2, V ,12, V ,12,' ,12) 

IF( NEVEN. NE. 1)  THEN 

C  WRITE  THE  ARRAY  PARAMETERS  USING  NON-EVEN  SYMMETRY. 

C  WRITE  THE  SHADING 

DO  100  JJ=1 , IM 
WRITE(2,503)  SHADE(JJ) 

503  FORMAT( F7. 5) 

100  CONTINUE 

C  WRITE  THE  ELEMENT  TOLERANCES. 

DO  110  LL=1,M 
IF(RER(LL).EQ.O)  THEN 
WRITE(2,507) 

507  FORMATC  0.00000,  0.00000') 

GOTO  110 
ENDIF 

DO  110  K.R=  1 ,  MM 

IF(U(KK,1).NE. RER(LL) )  GOTO  110 

WRITE( 2 , 504)  100.*DREAL(ELEMNT(KK) ) , DIMAG (ELEMNT( KK) ) 

504  FORMATC F8. 5, ' , ' ,F9. 5) 

110  CONTINUE 

C  WRITE  THE  COORDINATES. 

DO  130  JJ=1 , IM 
WRITE( 2 , 505)  X(JJ).YCJJ) 

505  FORMATC F5. 3 , ' , ' ,F5. 3) 

130  CONTINUE 
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WRITE(2,506)  WAVE , DOLAM 
506  FORMAT(F8.4/F8.4/'0' ) 

C  WRITE  THE  ARRAY  PARAMETERS  USING  EVEN  SYMMETRY. 

ELSE 

DO  140  J J=1 , IM+1 
140  WRITE(2,503)  SHADE(JJ) 

DO  150  LL= 1 , M 
IF(RER(LL) .EQ.O)  THEN 
WRITE(2,507) 

GOTO  150 
ENDIF 

DO  150  KK=1 ,MM 

IF(U(KK, 1 ) . NE. RER(LL) )  GOTO  150 

WRITE(2,504)  100.*  DRE AL ( ELEMNT ( KK ) ) , D IMAG  ( ELEMNT ( KK) ) 

150  CONTINUE 

DO  170  JJ=1 , IM+1 
WRITE(2,505)  X(JJ),Y(JJ) 

170  CONTINUE 

WRITE(2,506)  WAVE, DOLAM 
ENDIF 
RETURN 
END 

SUBROUTINE  OPTIMO(II) 

C  SUBROUTINE  TO  SELECT  FOUR  ELEMENTS  AT  A  TIME  FOR  AN  ARRAY  POSITION 
C  BASED  ON  THE  VALUES  OF  THE  KENDIG  TERMS  SPECIFIED. 

INTEGE R  Q,P,O(8,40) , TERM( 3 ) , U( 200 , 2 ) , RER( 200) , RERR( 200) , OSS( 20) 
INTEGER  T(4),BR(200),S 
DIMENSION  IMR( 200) 

REAL  K( 8)  , A(  2)  , V< 2)  , D< 40) , AMP( 200) , PHASE( 200) 

COMMON/ COMMOl/  M,MM,IM,IM1 

COMMON/COMM02/  U 

COMMON/ COMM03/  AMP, PHASE 

COMMON/ COMM3A/  BR 

COMMON/ C0MM04/  Q 

COMMON/ C0MM06/  OSS, TERM 

COMMON / COMM09 /  T 

COMMON/ C0MM10/  D 

COMMON/ C0MM12/  RER,RERR, IMR 

P=0 

S=1 

D1=0.0 

C  GENERATE  ALL  POSSIBLE  PERMUTATIONS  OF  4  ELEMENTS  FROM  THE 
C  ELEMENTS  REMAINING  IN  THE  POOL  AND  FIND  THE  40  BEST  PERMUTATIONS. 
DO  10  Kl=l,Q 
DO  20  Ll=l,Q 
IF  (L1.EQ.K1)  GOTO  20 
DO  30  Ml= 1 , Q 

IF  (M1.EQ.L1.0R.M1.EQ.K1)  GOTO  30 
DO  40  N1=1,Q 

IF  (N1.EQ.M1.0R.N1.EQ.L1.0R.N1.EQ.K1)  GOTO  40 
P=l+P 


GOTO  80 
ENDIF 

CALL  DMAX(D1 ,8) 

IF  (U1.LE.D1)  THEN 
DO  100  KK=1 ,8 
IF  (D(KK).EQ.Dl)  THEN 
D(KK)=U1 
0( 5,KK)=0( 1 , JJ) 

0(6,KK)=0(2, JJ) 

0( 7 ,KK)=0( 3 , JJ) 

0(8,KK)=0(4, JJ) 

GOTO  80 
ENDIF 

100  CONTINUE 
ENDIF 

80  CONTINUE 
IN=8 

C  BASED  ON  THE  VALUES  OF  THE  FIRST  TERM,  (OR  FIRST  TWO  TERMS), 

C  SELECT  FOUR  ELEMENTS  ON  THE  BASIS  OF  THE  LAST  TERM. 

200  W=1 . 0E12 
S=l+S 

DO  110  JJ=1 , IN 

CALL  SETU( S,U1,0(5,JJ),0(7,JJ),0(6,JJ),0(8,JJ)) 

IF  (Ul.LT.W)  THEN 
W=U1 

T( 1 )=0( 5 , JJ) 

T(2)=0(6, JJ) 

T(3)=0(7 ,JJ) 

T(4)=0(8, JJ) 

ENDIF 
1 10  CONTINUE 

C  FOR  THE  FIRST  SELECTION,  SPECIFY  THE  REFERENCE  POINT. 

IF  (II. EQ. 1)  CALL  SETK 

C  RESTORE  ALL  THOSE  ELEMENTS  WHICH  WERE  TEMPORARILY  REMOVED. 

CALL  PREPNSET 

C  STORE  BOTH  THE  SELECTED  ELEMENTS'  SEQUENTIAL  AND  ASSIGNED  NUMBERS 
C  IN  ARRAY  ORDER. 

RER(OSS( II) )=BR(T( 1 ) ) 

RER( OSS( II)+IM)=BR(T( 2 ) ) 

RER(0SS( II)+2*IM)=BR(T( 3) ) 

RER(OSS( II)+3*IM)=BR(T( 4) ) 

NN=0 

DO  120  LL=1 ,4 
DO  130  KK= 1 , MM 

IF(U(KK, 1 ) . EQ. BR(T( LL) ) )  THEN 
RERR(OSS(II)+(LL-l)*IM)=ttK 
IF(NN.EQ. 3)  GOTO  140 
NN=NN+1 
GOTO  120 
ENDIF 

130  CONTINUE 
120  CONTINUE 
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140  RETURN 
END 

SUBROUTINE  PERMSERIN 

C  ROUTINE  TO  WRITE  THE  ELEMENT  AND  ARRAY  PARAMETERS  TO  A  FILE 
C  TO  BE  READ  BY  PERMSER. 

INTEGER  RER( 200) ,RERR( 200) , IMR( 200) ,U( 200,2) 

REAL  AMP 1(200, 40) .PHASE  1(200, 40) ,X(70),Y(70) ,SHADE(70) 
REAL  THETA(40) 

COMMON/ C0MM01/  M,MM,IM,IM1 
C0MM0N/C0MM02/  U 
COMMON/ C0MM1 1/  AMP1 , PHASE 1 
COMMON/ COMM1 2/  RER.RERR, IMR 
COMMON/COMM13/  NEVEN,N1L ,N2L,LSEED 
COMMON/COMM 14/  WAVE , DO LAM 
COMMON/COMM15/  X,Y, SHADE 
WRITE(7,500)  DOLAM.WAVE ,M,MM,NEVEN, IM 

500  F0RMAT(F8.4,F8.4/I3, ' , 10, 1 ,0 , L ’ /I3/I2 , ’ , ' , 13) 

N=0 

IF(NEVEN.EQ. 1)  N=1 
DO  5  JJ=1 , IM+N 

5  WRITE(7,501)  X(JJ),Y(JJ) ,SHADE( JJ) 

501  FORMAT(F5. 3 , ' , ' ,F5. 3 , ' , ' ,F7. 5) 

DO  20  JJ=1 ,MM 

WRITE( 7 , 509)  U(JJ,1) 

509  FORMAT( 13) 

DO  20  LL= 1 , IM+N 

WRITE ( 7,504)  1 00. *AMP1 ( J J , LL) , PHASE  1 ( J J , LL) 

504  FORMAT(F8. 5 , ' , ' ,F9. 5) 

20  CONTINUE 

CALL  SET IMR 
DO  30  KK=1 ,MM 

WRITE(7 ,508)  IMR(KK) ,RERR(KK) 

508  FORMAT(I3, ' , ' ,13) 

30  CONTINUE 

READ( 3,600)  PHI 1, DELPHI, DELTHE 

600  FORMAT(3F10.5) 

WRITE( 7 , 505)  PHI 1 .DELPHI, DELTHE 

505  FORMAT( '0'/F6.2, ' , ' .F6.2/F6.2) 

TEM=( 360.0-PHI1 )/DELPHI 
NPHI=TEM 

TEM1=INT(TEM) / 2 

IF(TEM/2.0.NE.TEM1)  NPHI=TEM+1.0 
C  READ  IN  THE  THETA  VALUES  CONSTRAINED  FOR  EACH  PHI. 

DO  10  11=1, NPHI 
10  READ( 3,601)  THETA(II) 

601  FORMAT(F10. 5) 

DO  15  11=1, NPHI 

15  WRITE( 7 , 506)  THETA(II) 

506  FORMAT(F6.2) 

READ( 3,602)  ZLOLIM 

602  FORMAT(F12. 5) 


Wv 


v  V,' 

■vy^ 

■w-n 


v\%  .v 
‘-•V-V/ 


m 

S' 

•M- 


WRITE( 7 , 507 )  ABS(LSEED/ 10) .ZLOLIM 
507  F0RMAT(I9, ' , ' ,F8.2) 

RETURN 

END 

SUBROUTINE  PREPNSET 

C  ROUTINE  TO  RESTORE  THOSE  ELEMENTS  THAT  WERE  TEMPORARILY 
C  DELETED  FROM  THE  SELECTION  POOL  FOR  USE  IN  SELECTION  FOR  THE 
C  NEXT  POSITION. 

INTEGER  U( 200, 2) ,BR( 200) ,T( 4) 

COMMON/ COMMOl/  M,MM,IM,IM1 
COMMON/ COMM02/  U 
COMMON/ COMM3A/  BR 
COMMON/COMM09/  T 

C  PERMANENTLY  DELETE  THOSE  ELEMENTS  JUST  SELECTED. 

DO  10  JJ=1 ,MM 

IF  (U(JJ,1).EQ. BR(T( l)).OR.U(JJ,l).EQ. BR(T( 2) ) . OR. U(JJ,1).EQ. BR(T( 
13)). OR. U( JJ, 1 ) . EQ. BR(T( 4) ) )  THEN 
U(JJ,2)  — 1 
ENDIF 

10  CONTINUE 

DO  20  JJ=1 ,MM 

IF  (U( JJ,2) .EQ.-2)  U( JJ,2)=0 
20  CONTINUE 
RETURN 
END 

SUBROUTINE  REMOVE(II) 

C  ROUTINE  TO  REMOVE  ELEMENTS  THAT  ARE  AWAY  FROM  THE  CURRENT  MEAN 
C  VALUES  FOR  ALL  BUT  THE  FIRST  POSITION,  THESE  MEANS  ARE  THOSE  OF  THE 

C  ELEMENTS  OCCUPYING  THE  FIRST  POSITION. 

INTEGER  U( 200, 2) ,Q,P 
COMPLEX* 16  ELEMNT( 200) 

REAL  K( 8) 

COMMON/ COMMOl/  M,MM,IM,IM1 
COMMON/ COMM02/  U 
COMMON/ COMM04/  Q 
COMMON/ COMM05/  K,E 
COMMON/COMM16/  ELEMNT 
DATA  K(8)/0.0/ 

C  K( 5-8)  ARE  THE  PARAMETERS  THAT  CONTROL  THE  REMOVAL  OF  ELEMENTS 
C  FROM  THE  SELECTION  POOL  BASED  ON  THEIR  DISTANCE  FROM  THE  REFERENCE 
C  MEAN,  WITH  A  SCALE  OF  THE  REFERENCE  VARIANCE. 

K(5)=.25 
K(6)=. 10 
K(7)=,05 
1  F=>25.0+K(  8) 

C  P  INDICATES  WHETHER  ELEMENTS  HAVE  BEEN  REMOVED  DURING  AN  ITERATION 
C  OF  THE  REMOVAL  CYCLE. 

5  P=0 


IF  (U(JJ,2). NE. O.OR. U( J J , 1 ) . EQ . 0)  GOTO  10 

T1=(DREAL(ELEMNT(JJ))-K(3))**2 

T2=(DIMAG(ELEMNT(JJ))-K(4))**2 

IF  (SQRT(T1+T2).LE.F*E)  GOTO  10 

U(JJ,2)=-2 

P=l+P 

10  CONTINUE 
CALL  SETQ 

C  REPLACE  IS  CALLED  IF  THEIR  ARE  TOO  FEW  ELEMENTS  REMAINING  IN  THE 
C  SELECTION  POOL. 

IF  (Q.LT. 10)  CALL  REPLACE(&1) 

IF  (Q.LE.14)  GOTO  50 

C  IF  THIS  IS  THE  SELECTION  FOR  THE  FIRST  POSITION,  THE  MEAN  AND 
C  VARIANCE  OF  THE  GROUP  MUST  BE  DETERMINED  AFTER  ELEMENTS  HAVE 
C  BEEN  REMOVED. 

IF  (II. NE. 1)  GOTO  40 
IF  (P.EQ.O)  GOTO  40 
CALL  STAT(II) 

GOTO  5 

40  IF  (F.GT.1.75)  F=F-K(5) 

IF  (F.LE. 1.75. AND. F.GT. 1.0)  F=F-K(6) 

IF  (F.LE. 1.0)  F=F-K( 7 ) 

GOTO  5 

50  WRITE( 4 , 500)  F 

500  FORMAT(5X, 'FINAL  F  =  ’.F6.2) 

WRITE( 4,501)  Q 

501  FORMAT(5X,' NUMBER  OF  ELEMENTS  REMAINING  =  ',13) 

CALL  SORT 

RETURN 

END 

SUBROUTINE  REPLACE(*) 

C  ROUTINE  TO  RESTORE  ALL  ELEMENTS  THAT  WERE  TEMPORARILY  DELETED  BECAUSE 
C  10  OR  FEWER  ELEMENTS  WERE  LEFT.  ADDITIONALLY,  THE  INCREMENTS  ARE 
C  DECREASED. 

REAL  K( 8) 

INTEGER  U( 200,2) 

COMMON/ COMMOl/  M,MM,IM,IM1 

COMMON/ C0MM02/  U 

COMMON/ COMMO 5/  K,E 

K(8)=K(8)+5.0 

K( 5)=. 1 

K(6)=.075 

K(7)=.025 

DO  10  KK= 1 , MM 

IF  (U(KK,2) .EQ.-2)  U(KK,2)=0 
10  CONTINUE 
RETURN  1 
END 
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SUBROUTINE  RESPONSE(L) 

C  ROUTINE  TO  DETERMINE  EACH  ELEMENTS  AMPLITUDE  AND  PHASE  RESPONSE  BASED 
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C  ON  ITS  EQUIVALENT  CIRCUIT  VALUES  AND  THE  RADIATION  LOADING  FOR 
C  EACH  ARRAY  POSITION. 

INTEGER  U(200,2) ,OSS(40) 

REAL  RA(40),XA(40) 

REAL* 16  R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11 ,R12,R13 
REAL* 1 6  S(200,9) ,PI 
COMPLEX* 16  ELEMNT( 200) 

REAL  AMP 1(200, 40) .PHASE  1(200, 40) 

COMMON/COMMOl/  M,MM,IM,IM1 
COMMON/ COMM02/  U 
COMMON/ COMMO 6/  OSS, TERM 
COMMON/ COMM08/  RA,XA,F,SIZE 
COMMON/ COMM8 A/  S 
COMMON/ COMM11/  AMP1 , PHASE1 
COMMON/ COMM1 6/  ELEMNT 
PI=4 . 0Q0*QATAN( 1 . 0Q0 ) 

DO  20  JJ=1 ,MM 

R1=S(JJ,2)/(S(JJ,5)-S(JJ,4)) 

R2=S(JJ,3)/(2.0Q3*PI*S(JJ,2)*R1) 

R3=l.0QO/(R2*(2.OQO*PI*S(JJ,2))**2.0QO) 

R4=(S(JJ,6)/1.0Q12)-R2 
R5=l . 0Q03/S( JJ ,3) 

R6=S(JJ,7) 

R7=2. OQO*PI*QEXT( F ) 

R8=(QEXT( SIZE) )**2.0Q0 
R9=R8**2. 0Q0 

R10=1.0Q0-(R7**2.0Q0)*R3*R4+(R4/R2)-R7*R4*QEXT(XA(OSS(L)))*R9 
R1 l=R7*R4*(R5+QEXT(RA(OSS(L) ) )*R9/(R6**2.0Q0) ) 

R12=R8/ (R6*QSQRT(R10**2.0Q0+R1 1**2. OQO) ) 

R1 3=QATAN2D(R1 1 ,R10) 

C  AMPl  AND  PHASE  1  ARE  USED  AS  INPUT  TO  THE  PERMSER  FILE. 

C  ELEMNT  IS  USED  DURING  THE  SELECTION  PROCESS. 

AMPl ( JJ,OSS(L) )=SNGL(R12) 

PHASE1(JJ,0SS(L) )=SNGL(R13) 

IF  (U( JJ , 2 ) . EQ.-l )  GOTO  20 
ELEMNT( J J ) =DCMPLX(R 1 2 , R 1 3 ) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  SETIMR 

C  ROUTINE  TO  RANDOMIZE  THE  PERMUTATION  ARRAY  FOR  USE  IN  PERMSER 
C  PROGRAM. 

C  THIS  ROUTINE  APPEARS  IN  THE  NEWSTART  SUBROUTINE  OF  PERMSER. 

INTEGER  RER( 200) , RERR( 200) , IMR( 200) ,RERT( 200) , IMRT( 200) ,RERRT( 200) 
INTEGER  U(200,2) ,T(200) 

DIMENSION  NN( 200) 

COMMON/COMMOl/  M,MM,IM,IM1 
COMMON/COMM02/  U 
COMMON/COMM  12/  RER, RERR, IMR 
COMMON/COMM1 3/  NEVEN , NIL , N2L , LSEED 
DO  10  KK= 1 ,MM 
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NN(KK)=0 

T(KK)=U(KK,2) 

10  CONTINUE 
LL=1 
ITER=0 

1  ITER=ITER+1 

Y»FLOAT(MM)*RAN(LSEED)+1.0 

NN1=INT(Y) 

IF(LL.NE.l)  GOTO  21 
NN(LL)=NN1 
LL=LL+1 
GOTO  1 

21  IF(ITER.GT. 10*MM)  GOTO  92 
DO  22  KK= 1 , LL 
IF(NN1.EQ.NN(KK) )  GOTO  1 

22  CONTINUE 
NN(LL)=NN1 
LL=LL+1 

IF(LL.GT.MM)  GOTO  94 
GOTO  1 

92  ITER=5*MM 
GOTO  1 

C  THOSE  PERMUTATION  ARRAY  SLOTS  THAT  REPRESENT  WORKBENCH  POSITIONS 

C  MUST  BE  FILLED  BY  ELEMENTS  THAT  HAVE  NOT  BEEN  SELECTED. 

94  DO  39  KL=4*I..rl  ,MM 
DO  35  JJ=1 ,MM 
IF  (T(JJ).EQ.-l)  GOTO  35 
RERR(KL)=JJ 
T(JJ)=-1 
GOTO  39 

35  CONTINUE 

39  CONTINUE 

DO  36  JJ=1,MM 
IF(RERR( JJ) .EQ.O)  THEN 
DO  37  KK=1 ,MM 
IF  (T(KK).EQ.-l)  GOTO  37 
RERR( JJ)=KK 
T(KK)=-1 
GOTO  36 

37  CONTINUE 
ENDIF 

36  CONTINUE 

C  PLACE  EACH  ARRAY  IN  THE  SAME  RANDOM  ORDER. 

DO  30  KK= 1 , MM 
IMRT( KK) = IMR( NN(KK) ) 

RERT(KK)=RER(NN(KK) ) 

RERRT(KK)=RERR(NN(KK) ) 

30  CONTINUE 

DO  40  KK=l,MM 
RER( KK) =RERT( KK) 

IMR( KK)=IMRT(KK) 

RERR(KK)=RERRT(KK) 
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40  CONTINUE 
RETURN 
END 

SUBROUTINE  SETK 

C  ROUTINE  TO  STORE  THE  MEAN  VALUES  AND  VARIANCES  FOR  THE  FIRST  FOUR 
C  ELEMENTS  SELECTED  FOR  USE  AS  A  REFERENCE  POINT. 

INTEGER  T(4) 

REAL  K(8),A(2),V(2) ,AMP( 200) ,PHASE( 200) 

COMMON/ C0MM03/  AMP, PHASE 
COMMON/COMM05/  K,E 
COMMON/COMM07/  A,V 
COMMON/ C0MM09/  T 

C  COMPUTE  THE  MEANS  AND  VARIANCES  FOR  THE  FOUR  SELECTED  ELEMENTS. 

A( 1 )=0.0 
A( 2)=0. 0 
V(1)=0.0 
V(2)=0. 0 
DO  10  JJ=1 ,4 

A( 1)=A( 1)+AMP(T( JJ) )/4.0 
A(2)=A(2)+PHASE(T( JJ) )/4.0 
10  CONTINUE 

DO  20  JJ=1 ,4 

V(1)=V( 1)+(AMP(T(JJ))-A( 1))**2.0 
V( 2)=V( 2)+(PHASE(T( JJ) )-A( 2) )**2.0 
20  CONTINUE 

V( 1)=SQRT(V( l)/3.0) 

V(2)=SQRT(V(2)/3.0) 

WRITE(4,500) 

500  FORMAT( 5X, 'THE  MEANS  AND  VARIANCES  OF  THE  FIRST  FOUR  CHOSEN  ARE*) 
WRITE( 4,501)  100.*A( 1),A(2),100.*V(1),V(2) 

501  FORMAT(5X, 'AMPLITUDE  MEAN:  ',F12.8,'  PHASE  MEAN:  ' ,F12.8/7X, ' VARI 
1ANCE:  ' ,F12. 8, 4X, 'VARIANCE:  '.F12.8) 

C  INITIALIZE  THE  REFERENCE  MEANS  AND  VARIANCES. 

K(1)*V(1) 

K(2)-V(2) 

K(3)=A( 1) 

K(4)-A(2) 

E=SQRT(K( 1 )**2. 0+K( 2)**2. 0) 

RETURN 

END 

SUBROUTINE  SETQ 

C  ROUTINE  TO  KEEP  TRACK  OF  HOW  MANY  ELEMENTS  ARE  AVAILABLE  IN  THE 
C  SELECTION  POOL. 

INTEGER  U( 200,2) ,Q 
COMMON/COMMOl/  M,MM,IM,IM1 
COMMON/COMM02/  U 
COMMON/ C0MM04/  Q 
Qs0 

DO  10  JJ= 1 , MM 

IF  ( U( J J , 2 ) . EQ. 0. AND. U(JJ.l).NE.-l)  Q=l+Q 
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10  CONTINUE 
RETURN 
END 

SUBROUTINE  SORT 

C  ROUTINE  TO  RANK  THE  REMAINING  ELEMENTS,  BEFORE  THE  OPTIMO 
C  SEARCH,  ACCORDING  TO  THEIR  DISTANCES  FROM  THE  REFERENCE 
C  MEAN  VALUES. 

INTEGER  U( 200 , 2 ) , BR( 200) , D , BRT( 200) , RER( 200) , RERR( 200) 

COMPLEX* 16  ELEMNT( 200) 

CHARACTER* 1  NAME(40) 

REAL  AMP(200) ,PHASE(200) ,A(2) ,V(2) ,C(20) ,K(8) ,AMPT( 200) 

REAL  PHASET(200) ,RA(40) ,XA(40) 

COMMON/ COMMOl/  M,MM,IM,IM1 
COMMON/ COMMO 2/  U 
COMMON/ COMMO 3/  AMP, PHASE 
COMMON/COMM3A/  BR 
COMMON/ COMMO 5/  K,E 
COMMON/COMM6A/  NAME 
COMMON/ C0MM08/  RA,XA,F,SIZE 
COMMON/ COMM1 2/  RER.RERR.IMR 
COMMON/ COMM1 3/  NEVEN, NIL, N2L, LSEED 
COMMON/ COMM1 6/  ELEMNT 
COMMON/ COMM17/  C 
DATA  KL/0/ 

KL=KL+1 

D=0 

DO  10  JJ=1 ,MM 

IF  (U(JJ,2).NE.O.OR.U(JJ,1).EQ.O)  GOTO  10 
D=l+D 

C(D)=SQRT( (DREAL(ELEMNT( JJ) )-K(3) )**2+(DIMAG(ELEMNT( JJ) )-K(4) )**2) 
AMP(D)=DREAL( ELEMNT ( J J) ) 

PHASE(D)=DIMAG(ELEMNT( JJ) ) 

BR(D)=U( JJ, 1) 

10  CONTINUE 

DO  11  JJ=1,D 
CALL  CMAX(Cl.D) 

DO  12  KK= 1 , D 
IF(C(KK).EQ.C1)  THEN 
I=KK 
GOTO  16 
ENDIF 

IF(Cl.EQ.O.O)  GOTO  17 
12  CONTINUE 

16  AMPT(D+1-JJ)=AMP(I) 

PHASET(D+1-JJ)=PHASE( I) 

BRT(D+1-JJ)  =  BR(  I) 

C(I)=0.0 

1 1  CONTINUE 

17  DO  18  JJ=1 , D 
AMP( JJ)=AMPT( JJ) 

PHASE( JJ)=PHASET( JJ) 


BR( JJ)=BRT( JJ) 

18  CONTINUE 

C  DETERMINE  THE  MEANS  AND  PHASES  OF  THE  ELEMENTS  REMAINING  IN  THE  POOL. 
A( 1 )=0. 0 
A(2)=0.0 
V(1)=0.0 
V(2)=0. 0 
DO  20  JJ=1 , D 
A( 1 )=AMP( JJ)+A( 1 ) 

A(2)=PHASE(JJ)+A(2) 

20  CONTINUE 

A(  1 )=A( 1 )/FLOAT( D) 

A( 2)=A( 2)/FLOAT(D) 

DO  30  JJ=1 , D 

V(  1 )=( AMP( JJ)-A( 1 ) )**2+V( 1) 

V(2)=(PHASE( JJ)-A(2) )**2+V(2) 

30  CONTINUE 

V( 1 )=SQRT( V( 1 )/FLOAT(D-l ) ) 

V ( 2)=SQRT( V( 2)/FL0AT( D-l  )  ) 

WRITE(4, 102)  100.*A( 1 ) ,A(2) 

102  FORMAT(5X, 'AVERAGES  OF  CHOSEN  ELEMENTS:  AMPLITUDE  =  ',F12.8,3X,' 
l PHASE  =  ' ,F12.8) 

WRITE( 4,103)  100. *V( 1) ,V(2) 

103  FORMAT(5X, ’VARIANCES  OF  CHOSEN  ELEMENTS:  AMPLITUDE  =  ',F12.8,3X,' 
l PHASE  =  ' ,F12.8) 

C  IF  AN  ODD  ARRAY  IS  BEING  SELECTED  FOR,  THE  CENTER  ELEMENT  IS  THE 
C  ONE  WHICH  IS  CLOSEST  TO  THE  MEAN  AND  VARIANCE  OF  THE  REMAINING 
C  ELEMENTS. 

IF(KL.EQ.l.AND.NEVEN.EQ.l)  THEN 
RER( 1 )=BR( 1 ) 

DO  40  JJ=1 ,MM 

IF(U( JJ, 1 ) .EQ.RER( 1 ) )  THEN 
U( JJ,2)=-1 
RERR( 1 )=JJ 
ELSE 

U( JJ,2)=0 
ENDIF 

40  CONTINUE 

C  THE  CENTER  ELEMENTS  MEAN  AND  THE  GROUPS  VARIANCE  ARE  USED  AS  THE 
C  REFERENCE  POINT. 

K(3)=AMP(1) 

K(4)=PHASE( 1 ) 

K(  1)=V< l) 

K(2)*V(2) 

E=SQRT(K( 1)**2+K(2)**2) 

WRITE( 4,104)  100.*AMP( 1 ) , PHASE( 1 ) , NAME( 1 ) ,BR( 1 ) 

104  FORMAT(5X,' CENTER  ELEMENT  :  AMPLITUDE  =  ',F12.8,3X,'  PHASE  =  ' ,F 
112.8/'+' , 91X.A1 , 19X, 13) 

WRITE( 4,105)  RA( 1 ) , XA( 1 ) 

105  FORMAT(5X, 'RADIATION  LOADING:  REAL  =  ',E12.4,5X,'  IMAGINARY  = 

1 '  ,E12.4//) 

ENDIF 
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RETURN 

END 

SUBROUTINE  STAT(II) 

C  ROUTINE  TO  DETERMINE  THE  MEAN  VALUES  AND  VARIANCES  FOR  THE  ELEMENTS 
C  REMAINING  IN  THE  SELECTION  POOL. 

C  ADDITIONALLY,  IT  SPECIFIES  THE  REFERENCE  VALUES  IN  THE  SEARCH  FOR 
C  THE  ELEMENTS  TO  OCCUPY  THE  FIRST  ARRAY  POSITION. 

INTEGER  U( 200,2) , Q 
COMPLEX* 16  ELEMNT( 200) 

REAL  K(8),A(2),V(2) 

COMMON/ COMMOl/  M,MM,IM,IM1 

COMMON/COMM02/  U 

COMMON/ C0MM04/  Q 

COMMON/ C0MM05/  K,E 

COMMON/ COMMO 7/  A,V 

COMMON/ COMM1 6/  ELEMNT 

A( 1)=0.0 

A( 2)=0. 0 

V(1)=0.0 

V(2)=0.0 

CALL  SETQ 

DO  10  JJ=1 ,MM 

IF  (U(JJ,2).NE.O.OR.U(JJ,1).EQ.O)  GOTO  10 
A( 1)=DREAL(ELEMNT( JJ) )/FLOAT(Q)+A( 1 ) 

A(  2 )=DIMAG( ELEMNT( J J) ) / FLOAT ( Q)+A( 2 ) 

10  CONTINUE 

DO  20  JJ=1,MM 

IF  (U(JJ,2).NE.O.OR.U(JJ,1).EQ.O)  GOTO  20 
V( 1 )=(DREAL( ELEMNT ( JJ) )-A( 1 ) )**2+V( 1 ) 

V(2)=(DIMAG(ELEMNT( JJ) )-A( 2) )**2+V( 2) 

20  CONTINUE 

V( 1 )=SQRT( V( 1 )/FLOAT(Q-l ) ) 

V ( 2)=SQRT( V( 2)/FLOAT(Q-l ) ) 

IF  (II.GT.l)  RETURN 
K(1)=V(1) 

K(2)=V(2) 

K(3)=A( 1) 

K(4)=A(2) 

E=SQRT(K( 1 )**2+K( 2)**2) 

RETURN 

END 

SUBROUTINE  SETU( S , U1 , K2 , L2 ,M2 , N2) 

C  ROUTINE  TO  DETERMINE  THE  VALUE  OF  THE  KENDIG  TERM  SPECIFIED 
C  FOR  EACH  PERMUTATION  FOUND  IN  OPTIMO. 

C  K  CORRESPONDS  TO  POSITION  1 ,  L  TO  POSITION  3,  M  TO  POSITION  2 
C  AND  N  TO  POSITION  4. 

INTEGER  TERM( 3),S,0SS(40) 

REAL  AMP( 200) , PHASE( 200) 

COMMON/ COMMO 3/  AMP, PHASE 
COMMON/COMM06/  OSS, TERM 
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C  KENDIG  TERM  3. 

IF  (TERM(S).EQ.3)  THEN 

Ui=ABS( AMP(L2)*PHASE( L2) -AMP (K2)*PHASE(K2 )+AMP(N2)* PHASE (N2)-AMP(M 
12)*PHASE(M2) ) 

C  KENDIG  TERM  4. 

ELSEIF  (TERM(S).EQ.4)  THEN 

U1=ABS(AMP(L2)*PHASE(L2)-AMP(K2)*PHASE( K2)-AMP(N2)*PHASE(N2)+AMP(M 
12)*PHASE(M2) ) 

C  KENDIG  TERM  8. 

ELSE 

U1=ABS( AMP(L2 )*PHASE( L2)+AMP(K2 )* PHASE ( K2)-AMP( N2)*PHASE( N2)-AMP(M 
12)*PHASE(M2) ) 

ENDIF 

RETURN 

END 

SUBROUTINE  WRITEIT(L) 

C  ROUTINE  TO  WRITE  OUT  THE  RESULTS  FOR  EACH  POSITION. 

C  ADDITIONALLY,  A  MAP  OF  THE  NUMBERS  ASSIGNED  TO  EACH  ELEMENT 
C  IS  PRINTED  SHOWING  WHERE  IN  THE  ARRAY  EACH  ELEMENT  IS  LOCATED. 

INTEGER  U( 200, 2) ,BR( 200) ,RER(200) , IMR( 200) , NN( 200),OSS(40),T(4) 
INTEGER  RERR( 200) 

REAL  RA( 40) , XA( 40) 

CHARACTER* 1  NAME (40) 

COMPLEX* 16  ELEMNT(200) 

CHARACTER*  3  RRER( 200 ) , BLANK 
DATA  BLANK/'  '/ 

REAL  A(2) ,V(2) 

COMMON/COMMOl/  M,MM,IM,IM1 

COMMON/ COMMO 2/  U 

COMMON/ C0MM3 A/  BR 

COMMON/ C0MM06/  OSS, TERM 

COMMON/ COMM6 A/  NAME 

COMMON/ COMM08/  RA,XA,F,SIZE 

COMMON/ C0MM09/  T 

COMMON/ COMM1 2/  RER,RERR, IMR 

COMMON/ COMM1 3/  NEVEN.N1L, N2L, LSEED 

COMMON/ COMM1 6/  ELEMNT 

IF  (L.EQ.O)  GOTO  200 

WRITE(4,501)  NAME(0SS(L)),BR(T(1)),BR(T(2)),BR(T(3)),BR(T(4)) 

501  FORMAT( '+' ,91X,A1 , 19X, 13 , 3X, 13 , 3X, 13 , 3X, 13) 

WRITE(4,510)  RA(0SS(L) ) ,XA(OSS(L) ) 

510  FORMAT(5X, 'RADIATION  LOADING:  REAL  =  ' ,E12. 4, 5X, ' IMAGINARY  =  ' ,E1 
12.4//) 

RETURN 

C  DETERMINE  THE  MEANS  AND  VARIANCES  OF  ALL  THE  ELEMENTS  IN  THE  ARRAY. 
200  A( 1 )=0.0 
A(2)*0.0 
V(1)=0.0 
V(2)=0.0 
N=0 

DO  10  JJ=1 ,MM 
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IF  (U(JJ,2).NE.-1)  GOTO  10 
A( 1 )=DREAL( ELEMNT( JJ) )+A( 1 ) 

A(2)=DIMAG(ELEMNT(JJ))+A(2) 

N=N+1 

10  CONTINUE 

A( 1 )=A( 1 )/FLOAT(N) 

A( 2)=A( 2)/FLOAT(N) 

DO  20  JJ=1 ,MM 

IF  (U(JJ,2).NE.-1)  GOTO  20 

V( 1 )=( DREAL( ELEMNT( J J) )-A( 1 ) )**2+V( 1 ) 

V(2)=(DIMAG(ELEMNT( JJ) )-A( 2) )**2+V( 2) 

20  CONTINUE 

V( 1 )=SQRT(V( 1 )/FLOAT(N-l ) ) 

V( 2)=SQRT( V( 2)/FLOAT(N-l ) ) 

WRITE(4,502)  100. *A( 1 ) ,A( 2) 

502  FORMAT(5X, 'ARRAY  AMPLITUDE  MEAN  =  ',F12.8,':  ARRAY  PHASE  MEAN  = 

1' ,F12.8) 

WRITE(4 , 503)  100. *V( 1) ,V(2) 

503  FORMAT(5X, 'AMPLITUDE  VARIANCE  =  ’ ,F12.8, ' :  PHASE  VARIANCE  = 

1' ,F12.8) 

C  THE  ARRAY  RRER  IS  USED  TO  PRINT  A  MAP  OF  THE  ASSIGNED  NUMBERS 
C  OF  EACH  SELECTED  ELEMENT  IN  THE  POSITION  THAT  IT  OCCUPIES  IN 
C  THE  ARRAY. 

DO  36  JJ=1 , MM 
WRITE(RRER( JJ) , 505)  BLANK 
36  IMR(JJ)=JJ 
N=0 

IF(NEVEN.EQ. 1 )  THEN 
N=1 

WRITE( RRER( 1 ) , 504)  RER(l) 

ENDIF 

Nl=l+N 

DO  40  LL=1+N,IM 
DO  50  KK=1+N, IM1 
IF  (OSS(KK) .EQ. IMR(LL) )  THEN 
WRITE(RRER(N1) ,504)  RER(LL) 

WRITE(RRER(N1+IM) , 504)  RER(LL+IM) 

WRITE(RRER(N1+2*IM) ,504)  RER(LL+2*IM) 

WRITE(RRER(N1+3*IM) ,504)  RER( LL+3*IM) 

504  FORMAT(I3) 

N1=N1+1 
GOTO  40 

ENDIF 

50  CONTINUE 

WRITE(RRER( N1 ) , 505)  BLANK 
WRITE( RR£R( Nl+IM) , 505)  BLANK 
WRITE(RRER(N1+2*IM) , 505)  BLANK 
WRITE(RRER(N1+3*IM) ,505)  BLANK 

505  FORMAT(A3) 

Nl-Nl+1 

40  CONTINUE 

C  THIS  ROUTINE  INITIALIZES  AN  ARRAY,  NN,  TO  BE  USED  TO  PRINT  A  MAP  OF 
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C  THE  ELEMENT  NUMBERS  REPRESENTING  THEIR  POSITION  IN  THE  ARRAY. 

C  THE  SAME  ROUTINE  IS  FOUND  IN  THE  KENDANAL  SUBROUTINE  OF  THE  NEARPRES 
C  PROGRAM. 

LL=1 

IF(NEVEN.NE. 1)  THEN 

C  THIS  BRANCH  INITIALIZES  NN  FOR  AN  EVEN  ARRAY. 

NL1=2*IM 

NL2=IM 

DO  110  KK=1 , N2L 
DO  120  JJ=1 , NIL 
I1=JJ-1 

NN(LL)=IMR(NL1-I1 ) 

LL=LL+1 

120  CONTINUE 

NL1=NL1-N1L 

DO  130  JJ=N1L, 1,-1 

I1=JJ-1 

NN( LL) =IMR( NL2-I 1 ) 

LL=LL+1 

1 30  CONTINUE 

NL2-NL2-N1L 
1 10  CONTINUE 

NL1=2*IM+1 
NL2=3*IM+1 
DO  140  KK=1,N2L 
DO  150  JJ=N1L,1,-1 
IWJ-1 

NN(LL)=IMR(NL1+I1 ) 

LL=LL+1 

150  CONTINUE 

NL1=NL1+N1L 
DO  160  JJ=1,N2L 
I1=JJ-1 

NN( LL) =IMR( NL2+I 1 ) 

LL=LL+1 

160  CONTINUE 

NL2=NL2+N1L 
140  CONTINUE 

ELSE 

C  THIS  BRANCH  INITIALIZES  NN  FOR  AN  ODD  ARRAY,  (  AN  ARRAY  WITH 
C  AN  ELEMENT  AT  THE  CENTER. ) 

LL=1 

NL1=2MM+1 

NL2=IM+2-NlL 

DO  117  KK=1 , NIL 

DO  127  JJ=1 ,N1L*N2L,N1L 

Il-JJ-1 

NN(LL)=IMR(NL1-I1 ) 

127  LL=LL+1 

NL1=NL1-1 

DO  137  J J= 1 , NIL 

I1=.JJ-1 
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NN(LL)=IMR(NL2+I1 ) 

137  LL=LL+1 

NL2=NL2-N1L 
117  CONTINUE 

NL1=3*IM-(N2L-1 )*N1L+1 
DO  147  JJ=1 , NIL 
I1=JJ-1 

NN(LL)=IMR(NL1-I1) 

147  LL=LL+1 

NL1=NL1+N1L 
NN(LL)=IMR( 1) 

LL=LL+1 

DO  157  JJ=1 , NIL 
I1=JJ-1 

NN(LL)=IMR(NL2+I1) 

1 57  LL=LL+1 

NL2=4*IM+2-NlL*N2L 
DO  167  KK=1 , NIL 
DO  177  JJ=1 ,N1L 
I1=JJ-1 

NN( LL)=IMR( NL1-I1 ) 

177  LL=LL+1 

NL1=NL1+N1L 

DO  187  JJ=1 ,N1L*N2L,N1L 
I1=JJ-1 

NN( LL)=IMR( NL2+I 1 ) 

187  LL=LL+1 

NL2=NL2+1 
167  CONTINUE 

ENDIF 

C  THIS  STATEMENT  DETERMINES  WHETHER  THE  MAP  IS  TO  BE  PRINTED 
C  ON  A  NEW  PAGE,  BASED  ON  THE  SIZE  OF  THE  ARRAY. 

IF(M.LE. 36)  THEN 
WRITE(4,506) 

506  FORMAT( ' O' / ’ O' ,35X, ' PLACEMENT  OF  ELEMENTS  IN  THE  ARRAY.' /’O') 
ELSE 

WRITE(4, 509) 

509  FORMAT (' 1', 3 5X,' PLACEMENT  OF  ELEMENTS  IN  THE  ARRAY.' /'O') 

ENDIF 

DO  197  KK=1 ,M,N1L+N2L 
WRITE(4,507) 

WRITE( 4,508)  ( RRER( NN(LL) ) , LL=KK , KK+N2L+N1L- 1 ) 

WRITE(4,507) 

507  FORMAT ( IX) 

508  FORMAT( 35X,12(A3,3X)) 

197  CONTINUE 

RETURN 

END 
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APPENDIX  F 

LISTING  OF  ARRAY  RESPONSE  PROGRAM 


(]** ****** *************** ******* *********************************** 

C  NEARPRES  * 

C  DIRECTIONAL  BEAM  PATTERN  MAPPING  PROGRAM  * 

C  WRITTEN  BY  KING  W.  WIEMANN  * 

C* ******************** ***********  ********************************* 

DIMENSION  PHI( 50) , IDEG( 50) , THETA ( 200) ,AMPL( 50) , PHASE( 50) 
DIMENSION  SAMP( 50) ,AAMP( 50) , NAMP( 50) 

CHARACTER* 1  JAMP( 50) , ICOLON, IBLANK 

COMMON/ COMM0I/  X( 150) ,Y( 150) ,ELEMNT( 150) ,SHADE( 150) ,AK 

COMMON/ C0MM02/  ITER, NREPET , NPRINT ,K 

COMMON/ C0MM0 3/  THETA 1 , THETA2 , DE LTHE , PH 1 1 , PH 1 2 , DE LPH I 

COMMON/ COMM04/  D1,D2,DELD 

COMMON/ COMMO 5/  R 

COMMON/ C0MM06/  D 

COMMON/ COMM07 /  M 

COMMON/ COMM1 6/  NSOURC 

COMPLEX  ELEMNT 

DATA  ICOLON/ ' : ' / , IBLANK/ '  ' / 

L  =  0 

NCOUNT  =  0 

C  READ  TASK  INDEX  K,  AND  CALL  THE  APPROPRIATE  SUBROUTINE. 

READ( 3 , 10)  NSOURC, K.NTYPE 

10  FORMAT( 313) 

GOTO( 1 ,2, 3, 4, 5)  K 

1  CALL  ROUTE1(L,&20) 

2  CALL  ROUTE2(L,&20) 

3  CALL  ROUTE3(L,&20) 

4  CALL  ROUTE4(L,&20) 

5  CALL  ROUTE5(L,&20) 

C  IF  THE  TASK  CHOSEN  IS  TO  BE  DONE  ITERATIVELY  FOR  NEW  PARAMETERS 

C  DATAIN  IS  CALLED  FOR  EACH  SUCCESSIVE  ITERATION. 

C  AFTER  EACH  MAP  IS  PRINTED,  CONTROL  RETURNS  TO  LINE  11. 

11  CALL  DATAIN(L) 

IF( NPRINT. EQ. 1)  CALL  PRINT 
C  INITIALIZING  THE  THETA  AND  PHI  VALUES. 

C  FOR  K  =  1,2,3  THETA  IS  HELD  CONSTANT  WHILE  PHI  CYCLES. 

C  FOR  K  =  4,5,  PHI  IS  HELD  CONSTANT  WHILE  THETA  CYCLES. 

C 

C  DETERMINE  THE  RESPONSE  ON  THE  MAIN  BEAM.  (THETA  =  0.0; PHI  =  0.0) 
C 

20  IF(NTYPE.NE. 1)  GOTO  12 

CALL  PRESS  (0.0,0.0,AMPL1 , PHASE  1 ,DNORM,D) 

SAMP1  =  AMPL1/DNORM 
IF  (ABS(SAMPl).GT. 1.0592E-5)  THEN 
AAMP1  =  20. 0*ALOG10(SAMP1 ) 

ELSE 

AAMP1  =  99.0 
ENDIF 
Z=AAMP1 
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C  SET  UP  THE  PHI  AND  THETA  VALUES  FOR  TASKS  1,  2,  OR  3. 

12  IF(K.LE. 3)  THEN 

II  =  (PHI2  -  PHI1 )/DELPHI  +1.0 
IF  (II. GT. 37)  II  =  37 
JJ  =  ( THETA2  -  THETA1 )/DELTHE  +1.0 
IF  (JJ.GT.181)  JJ  =  181 
DO  25  I  =  1,11 
PHI(I)  =  PHI1  +  ( 1-1 )*DELPHI 
IDEG(I)  =  PHI(I) 

DO  30  J  =  1,JJ 

THETA(J)  =  THETA 1  +  ( J-l )*DELTHE 
C  PRINT  THE  HORIZONTAL  SCALE  FOR  THE  MAP. 

WRITE(4,35) 

35  FORMAT( ' 1 ' ,40X, 'MAP  OF  ARRAY  OUTPUT'/) 

WRITE  (4,40)  ( IDEG( I), 1=1, II, 2) 

40  FORMAT(5X, 'THETA' ,55X, 'PHI  (DEGREES) '/ 1 IX, 18(14 , 2X) , 14) 
ELSE 

C  SET  UP  THE  THETA  AND  PHI  VALUES  FOR  TASKS  4  AND  5. 


25 

30 


50 


55 


45 

60 


II  =  (THETA2  -  THETA1 )/DELTHE  +1.0 
IFQI.GT.37)  II  =  37 
JJ  =  (PHI 2  -  PHI1)/DELPHI  +  1.0 
IF(JJ.GT. 181)  JJ  =  181 
DO  55  I  =  1,11 

THETA(I)  =  THETA1  +  (I-1)*DELTHE 
IDEG(I)  =  THETA(I) 

PRINT  THE  HORIZONTAL  SCALE  FOR  THE  MAP. 

WRITE(4,35) 

WRITE  (4,45)  (IDEG(I) ,1=1,11,2) 

FORMAT(5X, 'PHI' ,55X, 'THETA  (DEGREES) '/ 1 IX, 18(14, 2X) , 14) 

DO  60  J  =  1,JJ 
PHI(J)  =  PHI1  +  ( J-l )*DELPHI 
ENDIF 

C  IF  A  HEIGHT  SPECTRUM  IS  THE  TASK,  THETA  AND  PHI  WILL  CYCLE  WHILE 
C  THE  VALUE  OF  D  IS  HELD  CONSTANT. 

65  IF(K.NE.5)  GOTO  67 

66  D=D1 

67  DO  80  J  =  1,JJ 
DO  70  I  =  1,11 
IF(K. GT. 3 )  THEN 

C  DEPENDING  ON  K,  DIFFERENT  ANGLE  PARAMETERS  ARE  SENT  TO  PRESS. 
CALL  PRESS  (THETA( I) , PHI( J) , AMPL( I) , PHASE( I ) , DNORM, D) 

ELSE 

CALL  PRESS  (THETA( J) , PHI( I) , AMPL( I ) , PHASE( I) , DNORM, D) 

ENDIF 

SAMP(I)  =  AMPL( I)/DNORM 
IF  (AflS(SAMP(I) ).GT. 1.0592E-5)  THEN 
AAMP(I)  =  -20. 0*ALOG10( SAMP( I) ) 

ELSE 

AAMP(I)  =  99.0 
ENDIF 

SUBRACT  THE  SIDE  LOBE  RESPONSE  FROM  THE  MAIN  BEAM  RESPONSE. 

IF( NTYPE. NE. 1 )  GOTO  78 
A AMP( I )=Z+AAMP( I) 

IF( AAMP( I) . GT. 99.5)  AAMP(I)=99.0 
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PHASE(I)=PHASE1-PHASE(I) 

78  NAMP(I)  =  INT( AAMP( I)  +  SIGN  (0. 5 ,AAMP( I) )) 

JAMP(I)  =  IBLANK 

IF<PHASE(I).GE. 90.0. AND. PRASE(I).LE. 270.0)  JAMP(I)  =  ICOLON 
IF(PHASE(I). LE.-90.0)  JAMP(I)  =  ICOLON 
70  CONTINUE 

IF(K.GT. 3)  THEN 

WRITE(4,84)  PHI( J) , (NAMP( I) , JAMP( I) ,1=1 ,11) 

ELSE 

WRITE  (4,84)  THETA( J) , ( NAMP( I) , JAMP( I) ,1=1,11) 

END  IF 

84  FORMAT  ( 5X,F6. 2 , 2X, 37( 12 ,A1 ) ) 

80  CONTINUE 

IF(K. NE. 5 )  GOTO  99 
D  =  D1  +  DELD 
IF(D.GT.D2)  GOTO  99 
GOTO  67 

99  WRITE  (4,110) 

110  FORMAT  (' O’, 5X,' THETA  IS  THE  ANGLE  MEASURED  FROM  THE  Z-AXIS  (ANGLE 

1  OF  INCIDENCE) '/5X, 'PHI  IS  THE  ANGLE  MEASURED  FROM  THE  X-AXIS  (RO 
2LL  PLANE) '/5X,' ALL  NUMBERS  ARE  NEGATIVE  DECIBELS' / IX, '**NOTE**  LAB 
3ELLING  OF  THE  HORIZONTAL  SCALE  IS  FOR  ALTERNATE  COLUMNS  AND  IS  TRU 
4NCATED  TO  INTEGER  MODE'/) 

GOTO( 101,102,103,104,105)  K 

101  WRITE( 4 , 1 30)  K 

130  FORMAT(5X, 'K=' 1X,I2,2X, ’THIS  IS  A  MAP  OF  THE  RESPONSE  OF  THE  ARRA 
1Y  TO  A  SIGNAL  INCIDENT  AT  ANGLES  PHI  AND  THETA'//) 

GOTO  136 

102  WRITE( 4,131)  K,R 

131  F0RMAT(5X,'K='1X,I2,2X, 'THIS  IS  A  MAP  OF  THE  PRESSURE  FIELD  ON  A  H 
1EMISHPERICAL  SURFACE' / IX, ' CENTERED  AT  THE  CENTER  OF  THE  ARRAY  AT  A 

2  DISTANCE  OF  R=' , 1X,F8.4//) 

GOTO  136 

103  WRITE(4, 132)  K,D 

132  FORMAT( 5X, 'K=' ,1X, I2,2X, 'THIS  IS  THE  PRESSURE  FIELD  IN  A  PLANE  PAR 
1ALLEL  TO  THE  PLANE  OF  THE  ARRAY  AT  HEIGHT  D= ' , IX, F8. 4//) 

GOTO  136 

104  WRITE(4, 133)  K,D,PHI1 

133  FORMAT( 5X, ' K=' , IX, 12 , 2X, ' THIS  IS  THE  PRESSURE  FIELD  AT  HEIGHT  D=', 
11X,F8.4,2X, 'AND  ALONG  THE  ROLL  PLANE  PHI=' , 1X,F8. 4//) 

GOTO  136 

105  WRITE(4, 134)  K, PHI, DI.D2, DELD 

134  FORMAT(5X, 'K=' 1X,I2,2X, 'THIS  IS  THE  PRESSURE  FIELD  AT  ROLL  PLANE  P 
1HI=' , IX, F8.4/5X, 'STARTING  AT  HEIGHT  Dl= ', IX, F8. 4 , 2X, ' ENDING  AT  HE I 
2GHT  D2=' ,1X,F8.4,2X,'AT  AN  INCREMENT  OF  DELTA  D=’ , IX, F8. 4//) 

C 

C  IF  THE  TASK  IS  TO  BE  REPEATED,  THE  COUNTER  IS  INCREMENTED  AND  THE  NEW 
C  VALUES  READ  IN. 

136  NCOUNT=NCOUNT+l 

IF(NCOUNT.GE. ITER)  GOTO  112 

L=NREPET 

GOTO  11 

C  IF  THE  PROGRAM  IS  FINISHED,  A  NEW  SET  OF  ARRAY  PARAMETERS  MAY  BE 


non 


C  CHANGED  AND  THE  PROGRAM  RUN  AGAIN. 

C  HOWEVER,  ONLY  ONE  PARAMETER  MAY  BE  CHANGED  AT  THIS  TIME. 

112  READ( 3, 137)  NAGAIN 

137  FORMAT(Il) 

IF( NAGAIN. NE. 1)  GOTO  113 
READ(3, 138)  ITER,NREPET,NPRINT 

138  FORMAT( 313) 

NC0UNT=0 
L=NREPET 
GOTO  11 

113  STOP 
END 

SUBROUTINE  PRESS( THETAP , PHIP , AMPLP , PHASEP , DNORM , D) 

C******* ******************************************************* ******* 

WRITTEN  BY  MARK  SCHAFER  * 

EXPANDED  BY  KING  W.  WIEMANN  * 

******************************************************************** 

COMMON/COMMO 1/  X( 1 50) , Y( 1 50) ,ELEMNT( 150), SHADE( 1 50) , AK 
COMMON/ COMMO 2/  ITER, NREPET, NPRINT, K 
COMMON/COMMO 5/  R 
COMMON/ COMMO 7/  M 
COMMON/ C0MM08/  RK 
DIMENSION  ANGLE( 150), A( 150) 

REAL* 8  U , SUMR , SUMI , DCOS , DS IN , DAT AN 2 , DAMPL 
COMPLEX  ELEMNT 

C  CONSTANTS  TO  CONVERT  DEGREES  TO  RADIANS  AND  VICE  VERSA. 

DATA  D2RAD/ 1.7453292 5  E-02/ 

C  INITIALIZE  THE  SUMS. 

SUMI  =0.0 
SUMR  =  SUMI 

UCOS  =  SIN( THETAP* D2RAD) *COS( PHIP*D2RAD) 

USIN  =  SIN(THETAP*D2RAD) *SIN( PHIP*D2RAD) 

C  SUM  THE  RESPONSE  OVER  ALL  THE  ELEMENTS. 

C  A(J)  IS  THE  ACTUAL  SHADING  OF  EACH  ELEMENT.  IT  IS  THE  PRODUCT  OF 
C  THE  IDEAL  SHADING  FOR  THE  DESIRED  ELEMENT  (A  FUNCTION  OF  POSITION), 

C  AND  THE  ELEMENTS  AMPLITUDE  RESPONSE. 

DNORM  =  0 
DO  10  J  =  1,M 

A(J)  =  REAL(ELEMNT( J) )*SHADE(J) 

DNORM  =  DNORM  +  SHADE(J) 

ANGLE(J)  =  AIMAG(ELEMNT( J) ) 

S=X( J)*UCOS+Y( J)*USIN 
U=DBLE( AK*S+ANGLE( J ) *D2RAD) 

C  DEPENDING  ON  THE  TASK  CHOSEN,  THE  RESPONSE  IS  COMPUTED  DIFFERENTLY. 

C  FOR  THE  PRESSURE  FIELD,  THERE  IS  A  PHASE  FACTOR  INVOLVING  THE  DISTANCE 
C  R  FROM  THE  CENTER  OF  THE  ARRAY.  FOR  K=l,  Rl=l  (i.e.  DIVISION  BY 
C  UNITY). 

GOTO( 5,4)  K 
R  =  D/ COS(THETAP) 
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R1  =  R  -  S 
GOTO  6 
R1  =  1.0 

DAMPL  =  DBLE( A( J) ) 

SUMR  =  SUMR  +  DC0S(U)*DAMPL/R1 
SUMI  =  SUMI  +  DSIN(U)*DAMPL/R1 
)  CONTINUE 

TAKE  THE  MAGNITUDE  OF  THE  AMPLITUDE  AND  CONVERT  THE  PHASE  TO  DEGREES. 
AMPLP  =  SNGL( SUMR* SUMR  +  SUMI*SUMI) 

AMPLP  =  SQRT( AMPLP) 

PHASEP  =  SNGL(DATAN2(SUMI ,SUMR) )/D2RAD 
IF(K.NE.l)  PHASEP  =  PHASEP-AK*R/D2RAD 
RETURN 
END 

SUBROUTINE  R0UTE1(L,*) 

COMMON/ COMMO 2/  ITER, NREPET, NPRINT,K 
COMMON/COMM03/  THETA1 , THETA2 , DELTHE , PHI 1 . PHI2 , DELPHI 
COMMON/ COMM 14/  NPOOL.MM 
READ( 3,10)  THETA 1.THETA2, DELTHE 
READ( 3,10)  PHI 1.PHI2, DELPHI 
READ( 3,15)  NPOOL 
0  FORMAT( 3F8. 4) 

5  FORMAT( 15) 

IF( NPOOL. EQ. 1)  READ( 3,15)  MM 
CALL  DATAIN(L) 

IF(NPRINT.EQ.l)  CALL  PRINT 

RETURN  1 

END 

SUBROUTINE  ROUTE2(L,*) 

COMMON/ C0MM02/  ITER, NREPET, NPRINT , K 

COMMON/ COMMO 3/  THETA1 ,THETA2 , DELTHE , PHI1 , PHI2 , DELPHI 
COMMON/ COMMO 5/  R 
READ( 3,10)  THETA1,THETA2, DELTHE 
READ( 3,10)  PHI1.PHI2, DELPHI 
READ( 3,15)  R 
0  FORMAT( 3F8. 4 ) 

5  FORMAT(F8.4) 

CALL  DATAIN(L) 

IF(NPRINT.EQ.  1)  CALL  PRINT 

RETURN  1 

END 

SUBROUTINE  ROUTE3(L,*) 

COMMON/ COMMO 2/  ITER, NREPET, NPRINT, K 
COMMON/ COMMO 3/  THETA1 , THETA2 , DELTHE, PHI1 , PHI2, DELPHI 
COMMON/ COMMO 6/  D 
READ( 3,10)  THETA1.THETA2, DELTHE 
READ(3, 10)  PHI1.PHI2, DELPHI 
READ( 3,15)  D 
0  FORMAT/ 3F8. 4) 
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15  FORMAT(F8.4) 

CALL  DATAIM(L) 

IF(NPRINT.EQ. 1)  CALL  PRINT 

RETURN  1 

END 

SUBROUTINE  ROUTE4(L,*) 

C0MM0N/C0MM02/  ITER, NREPET, NPRINT , K 
COMMON/ COMMO 3/  THETA1 , THETA2 , DELTHE , PHI 1 , PHI2 , DELPHI 
COMMON/COMM06/  D 
READ( 3,10)  THETA1,THETA2, DELTHE 
READ(3, 15)  PHI 
READ( 3,15)  D 
10  FORMAT(3F8.4) 

15  F0RMAT(F8. 4) 

PHI1  =  -PHI 
PHI2  =  PHI 
DELPHI  =  2*PHI 
CALL  DATAIN(L) 

IF(NPRINT.EQ. 1 )  CALL  PRINT 

RETURN  1 

END 

SUBROUTINE  ROUTE5(L,*) 

COMMON/ COMMO 2/  ITER , NREPET , NPRINT , K 

COMMON/ COMM03/  THETAl ,THETA2 , DELTHE, PHI 1 , PHI2 , DELPHI 
COMMON/ C0MM04/  D1,D2,DELD 
READ( 3,10)  THETA 1.THETA2, DELTHE 
READ(3, 15)  PHI 
READ( 3,10)  D1,D2,DELD 
10  FORMAT( 3F8. 4 ) 

15  F0RMAT(F8. 4) 

PHI1  =  -PHI 
PHI2  =  PHI 
DELPHI  =  2*PH1 
CALL  DATAIN(L) 

IF( NPRINT. EQ. 1)  CALL  PRINT 

RETURN  1 

END 

SUBROUTINE  PRINT 

C  ROUTINE  TO  PRINT  OUT  THE  ARRAY  PARAMETERS  FOR  THE  CURRENT  MAP. 
COMMON/COMMOl/  X( 150) ,Y( 150) ,ELEMNT( 150) ,SHADE( 150) ,AK 
C OMMO N / COMMO 2 /  ITER, NREPET , NPRINT ,K 
COMMON/ COMMO 3/  THETAl , THETA2 , DELTHE , PHI 1 , PHI2 , DELPHI 
COMMON/ COMM04/  D1 , D2 , DELD 
COMMGN/COMM05/  R 
COMMON/ COMMO 6/  D 
COMMON/ C0MM07/  M 
COMMON/ COMMO 9/  DOLAM.WAVE 
COMMON/ C0MM1 5/  NEVEN, NROUTE , IM , NIL , N2L 
COMPLEX  ELEMNT 
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DIMENSION  NN( 150) , IMR( 150) ,AMP( 150) ,PHASE( 150) 

DATA  L/0/ 

L=L+1 

DO  1  KK=1,M 

AMP(KK)=REAL(ELEMNT( KK) ) 

1  PHASE(KK)=AIMAG(ELEMNT(KK) ) 

IF(NROUTE.EQ. 1)  THEN 
IF(L.GT.l)  GOTO  153 
DO  4  KK= 1 , M 
4  IMR(KK)=KK 

LL=1 

C  THIS  BLOCK  INITIALIZES  AN  ARRAY,  NN,  WHICH  ALLOWS  THE  PRINTING 
C  ELEMENT  TOLERANCES  IN  THE  OUTPUT  FILE  SO  THAT  IT  MATCHES  THEIR 
C  RESPECTIVE  PLACEMENTS  IN  THE  ARRAY. 

IF(NEVEN.NE. 1)  THEN 

C  THIS  BRANCH  INITIALIZES  NN  FOR  AN  EVEN  ARRAY. 

NL1=2*IM 

NL2=IM 

DO  10  KK=1 , N2L 
DO  20  JJ=1 , NIL 
II=JJ-1 

NN(LL)=IMR(NL1-II) 

LL=LL+1 

20  CONTINUE 

NL1=NL1-N1L 

DO  30  JJ=N1L, 1,-1 

II=JJ-1 

NN(LL)=IMR(NL2-1I) 

LL=LL+1 

30  CONTINUE 

NL2=NL2-N1L 
10  CONTINUE 

NL1=2*IM+1 
NL2=3*IM+1 
DO  40  KK=1 , N2L 
DO  50  JJ=N1L ,1,-1 
II=JJ-1 

NN(LL)=IMR(NL1+II) 

LL=LL+1 

50  CONTINUE 

NL1=NL1+N1L 
DO  60  JJ=1 , N2L 
II=JJ-1 

NN(LL)=IMR(NL2+II) 

LL=LL+i 

60  CONTINUE 

NL2=NL2+N1L 
40  CONTINUE 

ELSE 

C  THIS  BRANCH  INITIALIZES  NN  FOR  AN  ODD  ARRAY,  (  AN  ARRAY  WITH 
C  AN  ELEMENT  AT  THE  CENTER.) 

LL=1 


OF 
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NL1=2*IM+1 

NL2=IM+2-NlL 

DO  17  KK= 1 , N 1 L 

DO  27  JJ=1 ,N1L*N2L,N1L 

II=JJ-1 

NN(LL)=IMR(NL1-II) 

LL=LL+1 
NL1=NL1-1 
DO  37  J J=1 , NIL 
II=JJ-1 

NN(LL)=IMR(NL2+II) 

LL=LL+1 

NL2=NL2-N1L 

CONTINUE 

NL1=3*IM-(N2L-1 )*N1L+1 
DO  47  JJ=1 , NIL 
II=JJ-1 

NN(  LL)=IMR( NL1-II) 

LL=LL+1 
NL1=NL1+N1L 
NN(LL)=IMR( 1) 

LL=LL+1 

DO  57  JJ=1 , NIL 
II=JJ-1 

NN(  LL) =IMR( NL2+I I) 

LL=LL+1 

NL2=4*IM+2-NlL*N2L 
DO  67  KK=1 , NIL 
DO  77  JJ=1 , NIL 
II=JJ-1 

NN(LL)=IMR(NL1-II) 

LL=LL+1 

NL1=NL1+N1L 

DO  87  JJ=1 , N1L*N2L , NIL 
II=JJ-1 

NN( LL)=IMR( NL2+II) 

LL=LL+1 

NL2=NL2+1 

CONTINUE 

ENDIF 

WRITE(4, 106) 

FORMAT( ’ 1 ' ,4X, 'PLACEMENT  OF  ELEMENT  TOLERANCES  IN  THE  ARRAY.') 
11=1 

DO  97  KK=1 ,M,N1L+N2L 
WRITE(4, 108) 

WRITE(4, 109)  (AMP(NN(LL) ) ,LL=KK,KK+N2L+N1L-1 ) 

WRITE(4, 108) 

WRITE(4, 109)  ( PHASE(NN( LL) ) , LL=KK,KK+N2L+N1L-1 ) 

WRITE(4, 108) 

FORMAT( IX) 

FORMAT(5X, 12(F9.5,2X)) 

11=11+1 
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97  CONTINUE 

CALL  KENDANAL 
ELSE 

C  THIS  WILL  PRINT  OUT  THE  ARRAY  PARAMETERS  FOR  NON-SYMMETRIC  ARRAYS. 

118  WRITE(4,245) 

245  FORMATS  1 SHADING,  POSITION(X, Y) ,  AND  TOLERANCES  (AMPL.  AND  PHASE 
1)'/ IX, 'FOR  EACH  ELEMENT') 

DO  279  J=1,M 

WRITE(4,250)  SHADE(J),  X(J),  Y(J),  ELEMNT(J) 

250  F0RMAT(1X,5(F9.5,2X)) 

279  CONTINUE 

291  END IF 

WRITE(4, 12)  DOLAM.WAVE 

12  FORMAT( 5X, ' D  OVER  LAMDA  = ' ,F8. 4 , 2X, ' WAVELENGTH= ' ,F8. 4///) 

C  EACH  TASK  HAS  A  UNIQUE  PARAMETER  ASSOCIATED  WITH  IT.  FOR  EACH 
C  TASK,  THAT  PARAMETER  IS  PRINTED. 

IF(K.NE. 2)  GOTO  21 
WRITE( 4,15)  K,R 

15  FORMAT( IX, 'K  =  ', 12 , 2X, ' THE  RADIUS  IS:  ’,F8.4///) 

21  IF(K.NE. 3)  GOTO  31 

WRITE(4, 16)  K,D 

16  FORMAT( IX, 'K  =  ',I2,2X,'THE  DISTANCE  IS:  ’ ,F8. 4///) 

31  IF(K.NE.4)  GOTO  41 

WRITE( 4,25)  K, PHI2 , D 

25  FORMAT ( IX, ' K  =  ' ,I2,2X,'THE  ROLL  PLANE  IS:  ' ,F8. 4 , 2X, 'THE  DISTANCE 

1  IS  :  ' ,F8. 4) 

41  IF(K.NE. 5)  GOTO  51 

WRITE(4 ,35)  K , PHI2 , D1 , D2 , DELD 

35  F0RMAT(1X,'K  =  ',I2,2X,'THE  ROLL  PLANE  IS:  ’ ,F8. 4/ IX, ' THE  INITIAL 
1AND  FINAL  DISTANCES  ARE:  ' , 2(F8. 4 , 2X) , ' AND  THE  INCREMENT  IS:  ' ,F8. 
24) 

51  RETURN 
END 

SUBROUTINE  KENDANAL 

C  ROUTINE  WHICH  DETERMINES  THE  KENDIG  TERMS,  AND  PLACES  THEM  IN 
C  THEIR  RESECTIVE  POSITION  IN  THE  FIRST  QUADRANT. 

COMMON/COMMOl/  X( 150) , Y( 1 50) , ELEMNT( 1 50) , SHADE( 1 50) ,AK 

COMMON/ COMMO 7/  M 

COMMON/ COMMIO/  IMR 

COMMON/ COMM 14/  NPOOL.MM 

COMMON/ COMM 1 5/  NEVEN, NROUTE , IM , NIL , N2L 

COMMON/COMM 16/  NSOURC 

COMPLEX  ELEMNT 

DIMENSION  V1(40),V2(40),V3(40),V4(40),V5(40),V6(40),V7(40) 
DIMENSION  V8(40),AMP( l 50) , PHASE( 1 50) , NN( 40) , IMR( 150) 

DATA  D2RAD/ 1.74532925  E-02/ 

DO  10  11=1, M 

AMP( 1 1 ) =REAL( ELEMNT( 1 1 ) ) 

PHASE( II)=AIMAG(ELEMNT( II) )*D2RAD 
10  CONTINUE 
KK1=IM 
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KK2=2*IM 

KK3=3*IM 

N=0 

IF(NEVEN.EQ. L)  N=1 
DO  20  KK= 1 , IM+N 
K1=KK+KK1 
K2=KK+KK2 
K3=KK+KK3 

V l 1=AMP(KK)*PHASE(KK) 

V22=AMP(K1 )*PHASE(K1 ) 

V33=AMP(K2)*PHASE(K2) 

V44=AMP(K3)*PHASE(K3) 

V1(KK)=AMP(KK)+AMP(K1)+AMP(K2)+AMP(K3) 

V  2 ( KK) = AM  P( K 1 ) +AMP ( K3  ) -AMP ( KK) -AMP ( K2 ) 

V3(KK)=V33-V1 1+V44-V22 
V4(KK)=V33-V1 1-V44+V22 

V  5 ( KK) =AMP( KK) -AMP( K2 ) +AMP ( K1)-AMP(K3) 
V6(KK)=AMP(KK;-AMP(K2)+AMP(K3)-AMP(K1 ) 

V7(KK)=V1 1+V22+V33+V44 

V8(KK)=V1 1+V33-V22-V44 
20  CONTINUE 

C  THE  ARRAY  NN  CONTROLS  THE  PLACEMENT  OF  THE  KENDIG  TERMS  SIMILAR  TO 
C  THE  WAY  NN  WORKS  IN  THE  PRINT  SUBROUTINE. 

IF(NEVEN.NE. 1 )  THEN 
C  INTITIALIZE  NN  FOR  AN  EVEN  ARRAY. 

LL=1 

NL1=N2L 

NL2=1 

DO  40  KK=1 , N2L 
DO  30  II=NL1 , NL2 , -1 
JJ=II-1 
NN(LL)=IM-JJ 
LL=LL+1 

30  CONTINUE 

NL1=NL1+N2L 
NL2=NL2+N2L 
40  CONTINUE 

ELSE 

C  INITIALIZE  NN  FOR  AN  ODD  ARRAY. 

LL=  1 
NT=1 

DO  60  JJ=1 , N2L 

DO  70  II=NT+N2L-1 ,NT,-1 

NN(LL)=IM+2-II 

LL=LL+l 

70  CONTINUE 

NT=NT+N1L 
60  CONTINUE 

DO  80  11=1, NIL 
NN( 1+(II-1)*N2L)=0 
80  CONTINUE 

ENDIF 
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NL1=1 
NL2=N2L 
WRITE(4, 109) 

109  FORMAT('r,15X,'KENDIG  TERMS'/) 

DO  50  11=1, N2L 
WR1TE(4, 104) 

IF(NEVEN.EQ. 1 . AND. NN(NL1 ) . EQ.O)  THEN 

WRITE(4, 106)  (V1(NN(JJ)),V2(NN(JJ)),JJ=NL1+1,NL2) 

WRITE(4, 106)  (V3(NN(JJ)),V4(NN(JJ)),JJ=NL1+1,NL2) 

WRITE(4, 106)  (V5(NN(JJ)), V6( NN( J J) ) , JJ=NL1+1 , NL2 ) 

WRITE(4, 106)  ( V7 (NN( JJ) ) , V8(NN( J J) ) , JJ=NL1+1 , NL2 ) 

ELSE 

WRITE( 4, 105)  ( V 1(NN( JJ) ) , V2(NN( JJ) ) , JJ=NL1 , NL2) 

WRITE(4, 105)  (V3(NN(JJ)),V4(NN(JJ)),JJ=NL1,NL2) 

WRITE( 4,105)  (V5(NN( JJ) ) , V6(NN( JJ) ) , J J=NL1 , NL2 ) 

WRITE(4, 105)  (V7(NN(JJ)),V8(NN(JJ)),JJ=NL1,NL2) 

ENDIF 

WRITE(4, 104) 

NL1=NL1+N2L 
NL2=NL2+N2L 
50  CONTINUE 

104  FORMAT (IX) 

105  FORMAT( 15X,5(2X,F8.5,2X,F8.5,2X)) 

106  FORMAT(37X,4(2X,F8.5,2X,F8.5,2X)) 

RETURN 

END 

SUBROUTINE  DATAIN(L) 

C  SUBROUTINE  TO  INPUT  DATA  FOR  NEARPRES. 

COMMON/ COMMOl/  X( 150) , Y( 150) ,ELEMNT( 150) .SHADE ( 150) ,AK 

COMMON/ COMMO 2/  ITER,NREPET,NPRINT,K 

COMMON/ COMMO 7/  M 

COMMON/ COMM08/  RK 

COMMON/COMM09/  DOLAM.WAVE 

COMMON/ COMMIO/  IMR 

COMMON/ COMM 14/  NPOOL.KM 

COMMON/ COMM1 5/  NEVEN,NR0UTE,IM,N1L,N2L 

COMMON/ COMM1 6/  NSOURC 

COMPLEX  ELEMNT 

DIMENSION  IMR( 150) 

PI  =  4 . 0*DATAN( 1 . 0D00) 

C  THIS  PASSES  CONTROL  TO  THE  SYMMETRIC  LOOP  ONLY  IF  DATAIN  IS  CALLED 
C  FROM  THE  MAIN  PROGRAM  RATHER  THAN  ONE  OF  THE  ROUTE  SUBROUTINES. 
IF(NROUTE.EQ. 1)  GOTO  500 
GOTO  (1,2, 3, 4)  L 

READ(3,95)  ITER,NREPET,NROUTE,M,NPRINT 
95  FORMAT( 513) 

IJ=M 

IF( NSOURC. EQ. 1)  IJ=MM 
C  NROUTE  SPECIFIES  LOOP  1  OR  LOOP  2 
IF(NROUTE.EQ. 1)  GOTO  500 
C 
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C  LOOP  ONE;  A  NONSYMMETRIC  ARRAY 
C 

C  READ  IN  THE  SHADING  COEFFICIENTS 

1  DO  10  J=1 ,M 

READ( 3,105)  SHADE(J) 

105  FORMAT(F9. 6) 

10  CONTINUE 

IF  (L.GT.O)  GOTO  50 
C  READ  IN  THE  ELEMENT  TOLERANCES 

2  IF( NSOURC. EQ. 1 )  THEN 

DO  93  J=1,IJ 
READ(3, 94 )  IMR(J) 

94  FORMAT(I4) 

93  CONTINUE 

ELSE 

DO  96  JJ=1,M 
IMR( JJ)=JJ 
96  CONTINUE 

ENDIF 

DO  15  JJ=1,IJ 

READ  ( 3, 1 iu)  ELEMNT(IMR(  JJ)) 

110  FORMAT (2F12.6) 

15  CONTINUE 

IF  (L.GT.O)  GOTO  50 
C  READ  IN  THE  ELEMENT  POSITIONS 

3  DO  20  J=1,M 

READ  (3,115)  X(J),  Y(J) 

115  FORMAT( 2F7 .3) 

20  CONTINUE 

IF  (L.GT.O)  GOTO  50 

C  READ  IN  THE  WAVELENGTH  AND  CONVERT  IT  TO  THE  WAVE  NUMBER  AK 

4  READ  (3,120)  WAVE 
120  FORMAT(F9.4) 

AK=2*PI/WAVE 
DOLAM=X( 2 ) -X( 1 ) /WAVE 
50  RETURN 
C 

LOOP  TWO;  A  SYMMETRIC  ARRAY. 

500  GOTO(401 ,402,403,404)  L 

READ( 3,405)  IM , NEVEN , N 1L , N2L 
405  FORMAT( 415) 

C  INPUT  THE  SHADING  COEFFICIENTS. 

C  IF  THE  ARRAY  IS  'ODD',  NEVEN=1 .  (’ODD’  ARRAY  MEANS  THERE  IS  A 
C  CENTER  ELEMENT,  AT  X=0,Y=0). 

401  IF(NEVEN.EQ. 1 )  THEN 

READ(3,415)  SHADE(l) 

DO  411  J = 2 , IM+ 1 
READ( 3,4 15)  SHADE(J) 

C  ASSIGN  EACH  POSITION  IN  THE  ARRAY  WITH  ITS  PROPER  SHADING. 
SHADE(J+IM)=SHADE(J) 

S HADE ( J + 2 * IM ) = S HADE ( J ) 
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SHADE( J+3*IM)=SHADE( J) 

411  CONTINUE 
ELSE 

C  INPUT  THE  VALUES  FOR  AN  EVEN  ARRAY. 

406  DO  410  J=1 , IM 

READ(3,415)  SHADE(J) 

C  ASSIGN  EACH  POSITION  IN  THE  ARRAY  WITH  ITS  PROPER  SHADING. 
SHADE(J+IM)=SHADE(J) 

SHADE (J+( 2*IM) )=SHADE( J) 

SHADE( J+( 3*IM) )=SHADE( J) 

415  FORMAT(F9.6) 

410  CONTINUE 

ENDIF 

IF  (L.GT.O)  GOTO  600 
C  READ  IN  THE  ELEMENT  TOLERANCES. 

402  IF( NSOURC. EQ. 1)  THEN 

DO  493  JJ=1,IJ 
READ(3,494)  IMR(JJ) 

494  FORMAT(I4) 

493  CONTINUE 

ELSE 

492  DO  496  JJ=1,IJ 

IMR(JJ)=JJ 
496  CONTINUE 

ENDIF 

DO  420  JJ=1,IJ 

READ  (3,110)  ELEMNT(IMR(JJ)) 

420  CONTINUE 

IF  (L.GT.O)  GOTO  600 
C  READ  IN  THE  ELEMENT  POSITIONS. 

C  IF  THE  ARRAY  IS  ODD,  NEVEN=1 .  (ODD  ARRAY  MEANS  THERE  IS  A  CENTER 
C  ELEMENT) . 

403  TF(NEVEN.EQ. 1)  THEN 

READ( 3,433)  X(1),Y(1) 

DO  438  J=2 , IM+1 
READ(3,433)  X(J),Y(J) 

C  ASSIGN  EACH  QUADRANT  WITH  THE  PROPER  SIGN  VALUES  FOR  X  AND  Y. 

Y( J+IM)=X( J) 

X( J+IM)=-Y( J) 

X( J+2*IM)=-X( J) 

Y( J+2*IM)=-Y( J) 

X( J+3*IM)=Y( J) 

Y( J+3*IM)=-X( J) 

438  CONTINUE 

ELSE 

C  INPUT  THE  POSITIONS  FOR  AN  EVEN  ARRAY. 

DO  430  J= 1 , IM 
READ(3,433)  X(J),Y(J) 

C  ASSIGN  EACH  QUADRANT  WITH  THE  PROPER  SIGN  VALUES  FOR  X  AND  Y. 
X(J+IM)=-X(J) 

X( J+( 2*IM) )=-X( J) 

X(J+(3*IM))=X(J) 
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Y(J+IM)=Y(J) 

Y(J+(2*IM))=-Y(J) 

Y(J+(3*IM))=-Y(J) 

33  FORHAT(2F7.3) 

30  CONTINUE 

ENDIF 

38  IF  (L.GT.O)  GOTO  600 

READ  IN  THE  WAVELENGTH  WAVE  AND  SPECIFY  THE  WAVENUMBER  AK. 

34  READ(3, 435)  WAVE 
READ(3,435)  DOLAM 

35  F0RMAT(F10.6) 

AK=2*PI*DOLAM 

RK=2*PI/WAVE 

30  RETURN 
END 
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